Ewald Method Research Articles

An integral equation method for the advection-diffusion equation on time-dependent domains in the plane

Boundary integral methods are attractive for solving homogeneous elliptic partial differential equations on complex geometries, since they can offer accurate solutions with a computational cost that is linear or close to linear in the number of discretization points on the boundary of the domain. However, these numerical methods are not straightforward to apply to time-dependent equations, which often arise in science and engineering. We address this problem with an integral equation-based solver for the advection-diffusion equation on moving and deforming geometries in two space dimensions. In this method, an adaptive high-order accurate time-stepping scheme based on semi-implicit spectral deferred correction is applied. One time-step then involves solving a sequence of non-homogeneous modified Helmholtz equations, a method known as elliptic marching. Our solution methodology utilizes several recently developed methods, including special purpose quadrature, a function extension technique and a spectral Ewald method for the modified Helmholtz kernel. Special care is also taken to handle the time-dependent geometries. The numerical method is tested through several numerical examples to demonstrate robustness, flexibility and accuracy.

Non-Ewald methods for evaluating the electrostatic interactions of charge systems: similarity and difference.

In molecular simulations, it is essential to properly calculate the electrostatic interactions of particles in the physical system of interest. Here we consider a method called the non-Ewald method, which does not rely on the standard Ewald method with periodic boundary conditions, but instead relies on the cutoff-based techniques. We focus on the physicochemical and mathematical conceptual aspects of the method in order to gain a deeper understanding of the simulation methodology. In particular, we take into account the reaction field (RF) method, the isotropic periodic sum (IPS) method, and the zero-multipole summation method (ZMM). These cutoff-based methods are based on different physical ideas and are completely distinguishable in their underlying concepts. The RF and IPS methods are "additive" methods that incorporate information outside the cutoff region, via dielectric medium and isotropic boundary condition, respectively. In contrast, the ZMM is a "subtraction" method that tries to remove the artificial effects, generated near the boundary, from the cutoff sphere. Nonetheless, we find physical and/or mathematical similarities between these methods. In particular, the modified RF method can be derived by the principle of neutralization utilized in the ZMM, and we also found a direct relationship between IPS and ZMM.

Integral equation analysis of multiport H‐plane microwave circuits by using 2D rectangular cavity Green's functions accelerated by the Ewald Method

An integral equation technique for the analysis of multiport H-plane microwave circuits composed of an unlimited number of arbitrarily shaped metallic and/or dielectric elements is proposed. The structure under study is divided into different regions after the application of the surface equivalence principle. Following this approach, the access waveguide ports are modelled by means of parallel plate Green's functions, whereas the central region is characterised by 2D rectangular cavity Green's functions. The use of this kind of Green's functions represents a novelty for solving H-plane integral equation problems. For the first time, the Ewald method has been employed in order to accelerate the computation of the 2D rectangular cavity Green's functions and their derivatives, needed in the integral equation. In addition, convergence studies show the numerical convenience of using the Ewald method to achieve a fast and accurate evaluation of these Green's functions. Results show that the analysis can be very efficient if the structure is segmented in a proper way based on fitting the cavity resonator within the geometry of the structure under analysis. In fact, all the walls of the structure that are coincident with the auxiliary resonator need not to be discretised during the numerical solution. Finally, several simulation examples of practical inductive microwave circuits are presented and discussed, showing a good agreement and higher computational efficiency when compared to results provided by commercial full-wave software tools and measurements.

A random batch Ewald method for charged particles in the isothermal-isobaric ensemble.

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.

Cohesion energy simulation of inorganic layered alkaline-earth fluorohalides

Within the framework of an interionic potential model, electrostatic and repulsive energies, lattice self-potentials for distinct ions, and the Madelung constant were calculated for several technologically important layered alkaline-earth fluorohalide scintillators belonging to the matlockite family and crystallizing with the PbFCl-type structure. The Coulombic term was computed using the conventional Ewald method, where the formulas were adjusted to optimize the computer calculation, and the overlap repulsive term was computed by means of an empirical approach based on the compressible ionic theory. The dispersive contribution was quantified using well-known formulations, and the polarization contribution was determined purely by geometric considerations based on both the size and anisotropic coordination of the highly polarizable halogen atom. In general, the obtained results were found to be in close agreement with the available data, except for the contribution of short-range interactions in the lattice. A quantitative explanation has been proposed to elucidate the differences detected. It was pointed out that the structural stability of these lamellar structures can be understood in terms of the anisotropic coordination of halogen anions, especially the high coordination of metal cations combined with their sizes. Our calculations led to an accurate evaluation of the cohesive energy, which, to the best of our knowledge, has never been measured before. Finally, the corresponding results will be useful for a better understanding of the chemical bonds and structural behavior of PbFCl-type compounds at high pressures.

A numerical study of lattice dynamics of single BaFX (X = Cl, Br, and I) crystals

The confirmed fundamental and practical prospects of inorganic matlockite materials crystallizing with the PbFCl-type structure require knowledge of their lattice vibrations and interatomic forces, which contribute to their layered morphology. In this work, we studied the lattice dynamics of BaFX, where X = Cl, Br, and I, in the framework of the harmonic shell model by fitting its multiple parameters to evaluate both long-wavelength Raman and infrared vibrational frequencies. The dynamical matrix was used to calculate the eigenmodes and the Ewald method was implemented to calculate the corresponding long-range part. Furthermore, the atomic displacements assigned to the normal modes were analyzed. The results are discussed and compared with experimental and previous theoretical data. The adopted model has been found to provide a much better reproduction of the measured frequencies compared to earlier simulations and first-principles calculations. The computed short-range force constants matched well with the cohesion and lamellar structure. This significant success can be extended to other compounds in the series and several complex ionic crystals.

Optimization of the Ewald method for calculating Coulomb interactions in molecular simulations.

Practical implementations of the Ewald method used to compute Coulomb interactions in molecular dynamics simulations are hampered by the requirement to truncate its reciprocal space series. It is shown that this can be mitigated by representing the contributions from the neglected reciprocal lattice vector terms as a simple modification of the real space expression in which the real and reciprocal space series have slightly different charge spreading parameters. This procedure, called the α' method, enables significantly fewer reciprocal lattice vectors to be taken than is currently typical for Ewald, with negligible additional computational cost, which is validated on model systems representing different classes of charged system, a CsI crystal and melt, water, and a room temperature ionic liquid. A procedure for computing accurate energies and forces based on a periodic sampling of an additional number of reciprocal lattice vectors is also proposed and validated by the simulations. The convergence characteristics of expressions for the pressure based on the forces and the potential energy are compared, which is a useful assessment of the accuracy of the simulations in reproducing the Coulomb interaction. The techniques developed in this work can reduce significantly the total computer simulation times for medium sized charged systems, by factors of up to ∼5 for those in the classes studied here.

Software for Evaluating Long-Range Electrostatic Interactions Based on the Ewald Summation and Its Application to Electrochemical Energy Storage Materials.

Electrochemical characteristics such as open-circuit voltage and ionic conductivity of electrochemical energy storage materials are easily affected, typically negatively, by mobile ion/vacancy ordering. Ordered phases can be identified based on the lattice gas model and electrostatic energy screening. However, the evaluation of long-range electrostatic energy is not straightforward because of the conditional convergence. The Ewald method decomposes the electrostatic energy into a real space part and a reciprocal space part, achieving a fast convergence in each. Due to its high computational efficiency, Ewald-based techniques are widely used in analyzing characteristics of electrochemical energy storage materials. In this work, we present software not only integrating Ewald techniques for two-dimensional and three-dimensional periodic systems but also combining the Ewald method with the lattice matching algorithm and bond valence. It is aimed to become a useful tool for screening stable structures and interfaces and identifying the ionic transport channels of cation conductors.

Improved Random Batch Ewald Method in Molecular Dynamics Simulations.

The random batch Ewald (RBE) is an efficient and accurate method for molecular dynamics (MD) simulations of physical systems at the nano/microscale. The method shows great potential to solve the computational bottleneck of long-range interactions, motivating a necessity to accelerate short-range components of the nonbonded interactions for a further speedup of MD simulations. In this work, we present an improved RBE method for the nonbonding interactions by introducing the random batch idea to constructing neighbor lists for the treatment of both the short-range part of the Ewald splitting and the Lennard-Jones potential. The efficiency that the novel neighbor list algorithm owes to the stochastic minibatch strategy can significantly reduce the total number of neighbors. We obtain the error estimate and convergence by theoretical analysis and implement the improved RBE method in the LAMMPS package. Benchmark simulations are performed to demonstrate the accuracy and stability of the algorithm. Numerical tests on computer performance by conducting large-scaled MD simulations for systems including up to 0.1 billion water molecules are run on massive clusters with up to 50000 CPU cores, demonstrating the attractive features such as the high parallel scalability and memory-saving of the method in comparison to the existing methods.

Molecular dynamic simulation of carbon dioxide, methane, and nitrogen adsorption on Faujasite zeolite

Removing impurities such as carbon dioxide and nitrogen from natural gas is a technical challenge and one of the major concerns in natural gas treatment process. In this study, adsorption of CH4, N2, and CO2 on the Faujasite (FAU) zeolite has been studied using molecular dynamics simulation at temperatures of 293, 308, and 323 K and pressures up to 1 MPa. COMPASS force field was used to model the interactions between zeolite and guest molecules. Ewald and atom-based summation methods were used for the calculation of electrostatic and van der Waals forces, respectively. Simulated results were modeled using Langmuir, Freundlich, Toth, and Sips adsorption isotherms. Sips isotherm for CO2, and Toth isotherm for CH4 and N2 pure compounds showed the best performance. Heat of adsorption for CH4, CO2, and N2 were calculated to be −15.48, −24.1, and −13.31 kJ⋅mol−1, respectively. A comparative study showed that the simulation model was successful in predicting the overall trend of the adsorption with acceptable accuracy.

Quantifying electrostatics from (bio)molecular dynamics trajectories using g-elpot

Electrostatics plays a pivotal role in many physiological processes, defining ion selectivity of ion channels, driving formation of supramolecular complexes, and determining energetics of transmembrane transporters, to name a few. Atomistic molecular dynamics (MD) simulations can serve as a rich source of structural and dynamic information on such processes. However, tools for robust high-resolution quantification of biomolecular electrostatics from MD trajectories are still limited. Here we present g_elpot, a GROMACS-based tool for calculating electrostatic potentials from MD simulations using the Smooth Particle Mesh Ewald (SPME) method.

Superscalability of the random batch Ewald method.

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 scale, but poor scalability limits the size of simulated systems. Here, we use 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 108 particles using 10 000 central processing unit 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, a micro-phase separated electrolyte, and a 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 promising solution on scalability of computing the Coulomb interaction. It is particularly useful and cost-effective to simulate ultra-large systems, which is either impossible or very costly to conduct using existing algorithms, and thus will be beneficial to a broad range of problems at nano-/micro-scales.

Sarkas: A fast pure-python molecular dynamics suite for plasma physics

We present an open-source, performant, pure-python molecular dynamics (MD) suite for non-ideal plasmas. The code, Sarkas, aims to accelerate the research process by providing an MD code complete with pre- and post-processing tools. Sarkas offers the ease of use of Python while employing the Numba library to obtain execution speeds comparable to that of compiled languages. The available tools in Sarkas include graphical displays of the equilibration process through a Jupyter interface and the ability to compute quantities such as, radial distribution functions, autocorrelation functions, and Green-Kubo relations. Many force laws used to simulate plasmas are included in Sarkas, namely, pure Coulomb, Yukawa and Molière pair-potentials. Sarkas also contains quantum statistical potentials and fast Ewald methods are included where necessary. An object-oriented approach allows for easy modification of Sarkas, such as adding new time integrators, boundary conditions and force laws. Program summaryProgram Title: SarkasCPC Library link to program files:https://doi.org/10.17632/zwpr5mpwms.1Developer's repository link:https://github.com/murillo-group/sarkasLicensing provisions: MITProgramming language: PythonNature of problem: Molecular dynamics (MD) is an important tool for non-ideal plasma physics research. The wealth of MD codes available are not designed for plasma physics problems. The available codes are written in low-level languages and do not provide pre- and post-processing libraries. These are instead written by researchers in interpreted languages, forcing researchers to have a high level of computing background.Solution method: Development of an MD suite for plasma physics, complete with pre- and post-processing tools most commonly used in plasma physics. The suite is entirely written in Python for enhanced user-friendliness. The slow speed of Python is mitigated by using the Numba library which is a just-in-time compiler for Python.

A hierarchical matrix approach for computing hydrodynamic interactions

For simulations of large numbers of small, spherical particles in a Stokes flow, the long-range hydrodynamic interactions approximated by the Rotne–Prager–Yamakawa (RPY) kernel can be summed rapidly using, for example, the fast multipole method (FMM) or the particle-mesh Ewald (PME) method. In this paper, we develop new fast methods for computing these sums using the H2 hierarchical matrix representation, for open and for periodic boundary conditions. To the best of our knowledge, the method for infinite periodic sums using the H2 hierarchical matrix representation is the first such method developed. We also consider a more general RPY kernel that handles polydisperse particle radii, and show analytically and experimentally that the proxy surface method for efficiently constructing the H2 hierarchical matrix representation remains effective in this case. Numerical tests demonstrate the well-controlled accuracy of the H2 summation methods and their linear-scaling computation and storage cost. We find that the H2 matrix approach has lower cost for computing the summations compared to FMM and PME, but higher precomputation cost (required for each particle configuration). This precomputation cost can be amortized over several summations when computing Brownian displacements or forces in Brownian and Stokesian dynamics simulations with very large numbers of particles.

Momentum removal to obtain the position-dependent diffusion constant in constrained molecular dynamics simulation.

The position-dependent diffusion coefficient along with free energy profile are important parameters needed to study mass transport in heterogeneous systems such as biological and polymer membranes, and molecular dynamics (MD) calculation is a popular tool to obtain them. Among many methodologies, the Marrink-Berendsen (MB) method is often employed to calculate the position-dependent diffusion coefficient, in which the autocorrelation function of the force on a fixed molecule is related to the friction on the molecule. However, the diffusion coefficient is shown to be affected by the period of the removal of the center-of-mass velocity, , which is necessary when performing MD calculations using the Ewald method for Coulombic interaction. We have clarified theoretically in this study how this operation affects the diffusion coefficient calculated by the MB method, and the theoretical predictions are proven by MD calculations. Therefore, we succeeded in providing guidance on how to select an appropriate value in estimating the position-dependent diffusion coefficient by the MB method. This guideline is applicable also to the Woolf-Roux method.

Validating the CHARMM36m protein force field with LJ-PME reveals altered hydrogen bonding dynamics under elevated pressures

The pressure-temperature phase diagram is important to our understanding of the physics of biomolecules. Compared to studies on temperature effects, studies of the pressure dependence of protein dynamic are rather limited. Molecular dynamics (MD) simulations with fine-tuned force fields (FFs) offer a powerful tool to explore the influence of thermodynamic conditions on proteins. Here we evaluate the transferability of the CHARMM36m (C36m) protein force field at varied pressures compared with NMR data using ubiquitin as a model protein. The pressure dependences of J couplings for hydrogen bonds and order parameters for internal motion are in good agreement with experiment. We demonstrate that the C36m FF combined with the Lennard-Jones particle-mesh Ewald (LJ-PME) method is suitable for simulations in a wide range of temperature and pressure. As the ubiquitin remains stable up to 2500 bar, we identify the mobility and stability of different hydrogen bonds in response to pressure. Based on those results, C36m is expected to be applied to more proteins in the future to further investigate protein dynamics under elevated pressures.

Partial Order–Disorder Transformation of Interpenetrated Porous Coordination Polymers

Controlling interpenetration and flexibility behaviors is intriguing and fundamental for porous coordination polymers. We report exceptional interpenetration behaviors involving controllable partia...

A fast spectral method for electrostatics in doubly periodic slit channels.

We develop a fast method for computing the electrostatic energy and forces for a collection of charges in doubly periodic slabs with jumps in the dielectric permittivity at the slab boundaries. Our method achieves spectral accuracy by using Ewald splitting to replace the original Poisson equation for nearly singular sources with a smooth far-field Poisson equation, combined with a localized near-field correction. Unlike existing spectral Ewald methods, which make use of the Fourier transform in the aperiodic direction, we recast the problem as a two-point boundary value problem in the aperiodic direction for each transverse Fourier mode for which exact analytic boundary conditions are available. We solve each of these boundary value problems using a fast, well-conditioned Chebyshev method. In the presence of dielectric jumps, combining Ewald splitting with the classical method of images results in smoothed charge distributions, which overlap the dielectric boundaries themselves. We show how to preserve the spectral accuracy in this case through the use of a harmonic correction, which involves solving a simple Laplace equation with smooth boundary data. We implement our method on graphical processing units and combine our doubly periodic Poisson solver with Brownian dynamics to study the equilibrium structure of double layers in binary electrolytes confined by dielectric boundaries. Consistent with prior studies, we find strong charge depletion near the interfaces due to repulsive interactions with image charges, which points to the need for incorporating polarization effects in understanding confined electrolytes, both theoretically and computationally.

Shadow Hamiltonian in classical NVE molecular dynamics simulations involving Coulomb interactions.

Microcanonical ensemble (NVE) Molecular Dynamics (MD) computer simulations are performed with negligible energy drift for systems incorporating Coulomb interactions and complex constraint schemes. In principle, such systems can now be simulated in the NVE ensemble for millisecond time scales, with no requirement for system thermostatting. Numerical tools for assessing drift in MD simulations are outlined, and drift rates of 10-6 K/μs are demonstrated for molten salts, polar liquids, and room temperature ionic liquids. Such drift rates are six orders of magnitude smaller than those typically quoted in the literature. To achieve this, the standard Ewald method is slightly modified so the first four derivatives of the real space terms go smoothly to zero at the truncation distance, rc. New methods for determining standard Ewald errors and the new perturbation errors introduced by the smoothing procedure are developed and applied, these taking charge correlation effects explicitly into account. The shadow Hamiltonian, Es, is shown to be the strictly conserved quantity in these systems, and standard errors in the mean of one part in 1010 are routinely calculated. Expressions for the shadow Hamiltonian are improved over previous work by accounting for O(h4) terms, where h is the MD time step. These improvements are demonstrated by means of extreme out-of-equilibrium simulations. Using the new methodology, the very low diffusion coefficients of room temperature 1-hexyl-3-methyl-imidazolium chloride are determined from long NVE trajectories in which the equations of motion are known to be integrated correctly, with negligible drift.

G_elpot: A Tool for Quantifying Biomolecular Electrostatics from Molecular Dynamics Trajectories.

Electrostatic forces drive a wide variety of biomolecular processes by defining the energetics of the interaction between biomolecules and charged substances. Molecular dynamics (MD) simulations provide trajectories that contain ensembles of structural configurations sampled by biomolecules and their environment. Although this information can be used for high-resolution characterization of biomolecular electrostatics, it has not yet been possible to calculate electrostatic potentials from MD trajectories in a way allowing for quantitative connection to energetics. Here, we present g_elpot, a GROMACS-based tool that utilizes the smooth particle mesh Ewald method to quantify the electrostatics of biomolecules by calculating potential within water molecules that are explicitly present in biomolecular MD simulations. g_elpot can extract the global distribution of the electrostatic potential from MD trajectories and measure its time course in functionally important regions of a biomolecule. To demonstrate that g_elpot can be used to gain biophysical insights into various biomolecular processes, we applied the tool to MD trajectories of the P2X3 receptor, TMEM16 lipid scramblases, the secondary-active transporter GltPh, and DNA complexed with cationic polymers. Our results indicate that g_elpot is well suited for quantifying electrostatics in biomolecular systems to provide a deeper understanding of its role in biomolecular processes.