Ewald Method Research Articles

The Spectral Ewald method for singly periodic domains

We present a fast and spectrally accurate method for efficient computation of the three dimensional Coulomb potential with periodicity in one direction. The algorithm is FFT-based and uses the so-called Ewald decomposition, which is naturally most efficient for the triply periodic case. In this paper, we show how to extend the triply periodic Spectral Ewald method to the singly periodic case, such that the cost of computing the singly periodic potential is only marginally larger than the cost of computing the potential for the corresponding triply periodic system.In the Fourier space contribution of the Ewald decomposition, a Fourier series is obtained in the periodic direction with a Fourier integral over the non-periodic directions for each discrete wave number. We show that upsampling to resolve the integral is only needed for modes with small wave numbers. For the zero wave number, this Fourier integral has a singularity. For this mode, we effectively need to solve a free-space Poisson equation in two dimensions. A very recent idea by Vico et al. makes it possible to use FFTs to solve this problem, allowing us to unify the treatment of all modes. An adaptive 3D FFT can be established to apply different upsampling rates locally. The computational cost for other parts of the algorithm is essentially unchanged as compared to the triply periodic case, in total yielding only a small increase in both computational cost and memory usage for this singly periodic case.

Kaiser-Bessel basis for particle-mesh interpolation.

In this work, we introduce the Kaiser-Bessel interpolation basis for the particle-mesh interpolation in the fast Ewald method. A reliable a priori error estimate is developed to measure the accuracy of the force computation in correlated charge systems, and is shown to be effective in optimizing the shape parameter of the Kaiser-Bessel basis in terms of accuracy. By comparing the optimized Kaiser-Bessel basis with the traditional B-spline basis, we demonstrate that the former is more accurate than the latter in part of the working parameter space, say, a relatively small real-space cutoff, a relatively small reciprocal space mesh, and a relatively large truncation of basis. In some cases, the Kaiser-Bessel basis is found to be more than one order of magnitude more accurate.

Publisher’s Note: Interacting dyon ensemble and confinement by particle mesh Ewald method [Phys. Rev. D 95 , 076006 (2017)

The free energy of a static quark-antiquark pair is obtained in an interacting dyon ensemble near the deconfinement temperature. Comparing the results with the noninteracting case, we observe that the string tension between the quark-antiquark pair increases for the interacting ensemble. As a result, the confinement temperature decreases.

Interacting dyon ensemble and confinement by particle mesh Ewald method

The free energy of a static quark-antiquark pair is obtained in an interacting dyon ensemble near the deconfinement temperature. Comparing the results with the noninteracting case, we observe that the string tension between the quark-antiquark pair increases for the interacting ensemble. As a result, the confinement temperature decreases.

Fast Ewald summation for free-space Stokes potentials

We present a spectrally accurate method for the rapid evaluation of free-space Stokes potentials, i.e., sums involving a large number of free space Green’s functions. We consider sums involving stokeslets, stresslets and rotlets that appear in boundary integral methods and potential methods for solving Stokes equations. The method combines the framework of the Spectral Ewald method for periodic problems (Lindbo and Tornberg in J Comput Phys 229(23):8994–9010, 2010. doi:10.1016/j.jcp.2010.08.026), with a very recent approach to solving the free-space harmonic and biharmonic equations using fast Fourier transforms (FFTs) on a uniform grid (Vico et al. in J Comput Phys 323:191–203, 2016. doi:10.1016/j.jcp.2016.07.028). Convolution with a truncated Gaussian function is used to place point sources on a grid. With precomputation of a scalar grid quantity that does not depend on these sources, the amount of oversampling of the grids with Gaussians can be kept at a factor of two, the minimum for aperiodic convolutions by FFTs. The resulting algorithm has a computational complexity of O(N log N) for problems with N sources and targets. Comparison is made with a fast multipole method to show that the performance of the new method is competitive.

A generalized Ewald decomposition for screened Coulomb interactions.

Medium-range interactions occur in a wide range of systems, including charged-particle systems with varying screening lengths. We generalize the Ewald method to charged systems described by interactions involving an arbitrary dielectric response function ϵ(𝐤). We provide an error estimate and optimize the generalization to find the break-even parameters that separate a neighbor list-only algorithm from the particle-particle particle-mesh algorithm. We examine the implications of different choices of the screening length for the computational cost of computing the dynamic structure factor. We then use our new method in molecular dynamics simulations to compute the dynamic structure factor for a model plasma system and examine the wave-dispersion properties of this system.

Anisotropy in geometrically rough structure of ice prismatic plane interface during growth: Development of a modified six-site model of H2O and a molecular dynamics simulation.

This paper presents a modified version of the six-site model of H2O [H. Nada and J. P. J. M. van der Eerden, J. Chem. Phys. 118, 7401 (2003)]. Although the original six-site model was optimized by assuming the cut-off of the Coulomb interaction at an intermolecular distance of 10 Å, the modified model is optimized by using the Ewald method for estimating the Coulomb interaction. Molecular dynamics (MD) simulations of an ice-water interface suggest that the melting point of ice at 1 atm in the modified model is approximately 274.5 K, in good agreement with the real melting point of 273.15 K. MD simulations of bulk ice and water suggest that the modified model reproduces not only the structures and density curves of ice and water, but also the diffusion coefficient of water molecules in water near the melting point at 1 atm. Using the modified model, a large-scale MD simulation of the growth at an ice-water interface of the prismatic plane is performed to elucidate the anisotropy in the interface structure during growth. Simulation results indicate that the geometrical roughness of the ice growth front at the interface is greater in the c-axis direction than in the direction normal to the c-axis when it is analyzed along the axes parallel to the prismatic plane. In addition, during the growth at the interface, the transient appearance of specific crystallographic planes, such as a {202¯1} pyramidal plane, occurs preferentially at the ice growth front. The effect of different ensembles with different simulation systems on the anisotropy in the interface structure is also investigated.

1-D Periodic Green’s Function for Leaky and Complex Waves Using the Ewald Method

The Ewald method for evaluating the free-space 1-D periodic Green’s function is extended to leaky waves by allowing for complex wavenumbers. It is shown that care must be taken when choosing the path of integration in the complex plane in the derivation of the Ewald method in order to obtain solutions that correspond to physical leaky-wave solutions. The use of different paths results in an evaluation of the Ewald method using an analytical continuation of the exponential integral function previously used when the wavenumber is real. The extension of the Ewald method to complex wavenumbers allows for the treatment of practical periodic leaky-wave antennas and periodic guiding structures with losses, as well as metamaterials, including 1-D chains of plasmonic nanoparticles.

Multiple Staggered Mesh Ewald: Boosting the Accuracy of the Smooth Particle Mesh Ewald Method.

The smooth particle mesh Ewald (SPME) method is the standard method for computing the electrostatic interactions in the molecular simulations. In this work, we develop the multiple staggered mesh Ewald (MSME) method, which averages the SPME forces computed on, for example, M, staggered meshes. We prove, from a theoretical perspective, that the MSME is as accurate as the SPME, but uses M2 times less mesh points in a certain parameter range. In the complementary parameter range, the MSME is as accurate as the SPME with twice the interpolation order. The theoretical conclusions are numerically validated both by a uniform and uncorrelated charge system, and by a three-point-charge water system that is widely used as solvent for the biomacromolecules.

Isotropic periodic sum for multipole interactions and a vector relation for calculation of the Cartesian multipole tensor.

Isotropic periodic sum (IPS) is a method to calculate long-range interactions based on the homogeneity of simulation systems. By using the isotropic periodic images of a local region to represent remote structures, long-range interactions become a function of the local conformation. This function is called the IPS potential; it folds long-ranged interactions into a short-ranged potential and can be calculated as efficiently as a cutoff method. It has been demonstrated that the IPS method produces consistent simulation results, including free energies, as the particle mesh Ewald (PME) method. By introducing the multipole homogeneous background approximation, this work derives multipole IPS potentials, abbreviated as IPSMm, with m being the maximum order of multipole interactions. To efficiently calculate the multipole interactions in Cartesian space, we propose a vector relation that calculates a multipole tensor as a dot product of a radial potential vector and a directional vector. Using model systems with charges, dipoles, and/or quadrupoles, with and without polarizability, we demonstrate that multipole interactions of order m can be described accurately with the multipole IPS potential of order 2 or m - 1, whichever is higher. Through simulations with the multipole IPS potentials, we examined energetic, structural, and dynamic properties of the model systems and demonstrated that the multipole IPS potentials produce very similar results as PME with a local region radius (cutoff distance) as small as 6 Å.

Rapid QM/MM approach for biomolecular systems under periodic boundary conditions: Combination of the density-functional tight-binding theory and particle mesh Ewald method.

A quantum mechanical/molecular mechanical (QM/MM) approach based on the density-functional tight-binding (DFTB) theory is a useful tool for analyzing chemical reaction systems in detail. In this study, an efficient QM/MM method is developed by the combination of the DFTB/MM and particle mesh Ewald (PME) methods. Because the Fock matrix, which is required in the DFTB calculation, is analytically obtained by the PME method, the Coulomb energy is accurately and rapidly computed. For assessing the performance of this method, DFTB/MM calculations and molecular dynamics simulation are conducted for a system consisting of two amyloid-β(1-16) peptides and a zinc ion in explicit water under periodic boundary conditions. As compared with that of the conventional Ewald summation method, the computational cost of the Coulomb energy by utilizing the present approach is drastically reduced, i.e., 166.5 times faster. Furthermore, the deviation of the electronic energy is less than 10-6 Eh. © 2016 Wiley Periodicals, Inc.

Pressure Profile Calculation with Mesh Ewald Methods.

The importance of calculating pressure profiles across liquid interfaces is increasingly gaining recognition, and efficient methods for the calculation of long-range contributions are fundamental in addressing systems with a large number of charges. Here, we show how to compute the local pressure contribution for mesh-based Ewald methods, retaining the typical N log N scaling as a function of the lattice nodes N. This is a considerable improvement on existing methods, which include approximating the electrostatic contribution using a large cutoff and the, much slower, Ewald calculation. As an application, we calculate the contribution to the pressure profile across the water/vapor interface, coming from different molecular layers, both including and removing the effect of thermal capillary waves. We compare the total pressure profile with the one obtained using the cutoff approximation for the calculation of the stresses, showing that the stress distributions obtained using the Harasima and Irving-Kirkwood path are quite similar and shifted with respect to each other at most 0.05 nm.

Crystal structure prediction of rigid molecules.

A non-polarizable force field based on atomic multipoles fit to reproduce experimental crystal properties and ab initio gas-phase dimers is described. The Ewald method is used to calculate both long-range electrostatic and 1/r(6) dispersion energies of crystals. The dispersion energy of a crystal calculated by a cutoff method is shown to converge slowly to the exact Ewald result. A method for constraining space-group symmetry during unit-cell optimization is derived. Results for locally optimizing 4427 unit cells including volume, cell parameters, unit-cell r.m.s.d. and CPU timings are given for both flexible and rigid molecule optimization. An algorithm for randomly generating rigid molecule crystals is described. Using the correct experimentally determined space group, the average and maximum number of random crystals needed to find the correct experimental structure is given for 2440 rigid single component crystals. The force field energy rank of the correct experimental structure is presented for the same set of 2440 rigid single component crystals assuming the correct space group. A complete crystal prediction is performed for two rigid molecules by searching over the 32 most probable space groups.

Computationally efficient analysis of extraordinary optical transmission through infinite and truncated subwavelength hole arrays.

The authors present a computationally efficient technique for the analysis of extraordinary transmission through both infinite and truncated periodic arrays of slots in perfect conductor screens of negligible thickness. An integral equation is obtained for the tangential electric field in the slots both in the infinite case and in the truncated case. The unknown functions are expressed as linear combinations of known basis functions, and the unknown weight coefficients are determined by means of Galerkin's method. The coefficients of Galerkin's matrix are obtained in the spatial domain in terms of double finite integrals containing the Green's functions (which, in the infinite case, is efficiently computed by means of Ewald's method) times cross-correlations between both the basis functions and their divergences. The computation in the spatial domain is an efficient alternative to the direct computation in the spectral domain since this latter approach involves the determination of either slowly convergent double infinite summations (infinite case) or slowly convergent double infinite integrals (truncated case). The results obtained are validated by means of commercial software, and it is found that the integral equation technique presented in this paper is at least two orders of magnitude faster than commercial software for a similar accuracy. It is also shown that the phenomena related to periodicity such as extraordinary transmission and Wood's anomaly start to appear in the truncated case for arrays with more than 100 (10×10) slots.

Ambient-Potential Composite Ewald Method for ab Initio Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulation.

A new approach for performing Particle Mesh Ewald in ab initio quantum mechanical/molecular mechanical (QM/MM) simulations with extended atomic orbital basis sets is presented. The new approach, the Ambient-Potential Composite Ewald (CEw) method, does not perform the QM/MM interaction with Mulliken charges nor electrostatically fit charges. Instead the nuclei and electron density interact directly with the MM environment, but in a manner that avoids the use of dense Fourier transform grids. By performing the electrostatics with the underlying QM density, the CEw method avoids self-consistent field instabilities that have been encountered with simple charge mapping procedures. Potential of mean force (PMF) profiles of the p-nitrophenyl phosphate dissociation reaction in explicit solvent are computed from PBE0/6-31G* QM/MM molecular dynamics simulations with various electrostatic protocols. The CEw profiles are shown to be stable with respect to real-space Ewald cutoff, whereas the PMFs computed from truncated and switched electrostatics produce artifacts. PBE0/6-311G**, AM1/d-PhoT, and DFTB2 QM/MM simulations are performed to generate two-dimensional PMF profiles of the phosphoryl transesterification reactions with ethoxide and phenoxide leaving groups. The semiempirical models incorrectly produce a concerted ethoxide mechanism, whereas PBE0 correctly produces a stepwise mechanism. The ab initio reaction barriers agree more closely to experiment than the semiempirical models. The failure of Mulliken-charge QM/MM-Ewald is analyzed.

Extended Ewald summation technique

We present a technique to improve the accuracy and to reduce the computational labor in the calculation of long-range interactions in systems with periodic boundary conditions. We extend the well-known Ewald method by using a linear combination of screening Gaussian charge distributions instead of only one. This enables us to find faster converging real-space and reciprocal space summations. The combined simplicity and efficiency of our method is demonstrated, and the scheme is readily applicable to large-scale periodic simulations, classical as well as quantum. Moreover, apart from the required a priori optimization the method is straightforward to include in most routines based on the Ewald method within, e.g., density-functional, molecular dynamics, and quantum Monte Carlo calculations.

A critical appraisal of the zero-multipole method: Structural, thermodynamic, dielectric, and dynamical properties of a water system.

We performed extensive and strict tests for the reliability of the zero-multipole (summation) method (ZMM), which is a method for estimating the electrostatic interactions among charged particles in a classical physical system, by investigating a set of various physical quantities. This set covers a broad range of water properties, including the thermodynamic properties (pressure, excess chemical potential, constant volume/pressure heat capacity, isothermal compressibility, and thermal expansion coefficient), dielectric properties (dielectric constant and Kirkwood-G factor), dynamical properties (diffusion constant and viscosity), and the structural property (radial distribution function). We selected a bulk water system, the most important solvent, and applied the widely used TIP3P model to this test. In result, the ZMM works well for almost all cases, compared with the smooth particle mesh Ewald (SPME) method that was carefully optimized. In particular, at cut-off radius of 1.2 nm, the recommended choices of ZMM parameters for the TIP3P system are α ≤ 1 nm(-1) for the splitting parameter and l = 2 or l = 3 for the order of the multipole moment. We discussed the origin of the deviations of the ZMM and found that they are intimately related to the deviations of the equilibrated densities between the ZMM and SPME, while the magnitude of the density deviations is very small.

Multilevel summation with B-spline interpolation for pairwise interactions in molecular dynamics simulations.

The multilevel summation method for calculating electrostatic interactions in molecular dynamics simulations constructs an approximation to a pairwise interaction kernel and its gradient, which can be evaluated at a cost that scales linearly with the number of atoms. The method smoothly splits the kernel into a sum of partial kernels of increasing range and decreasing variability with the longer-range parts interpolated from grids of increasing coarseness. Multilevel summation is especially appropriate in the context of dynamics and minimization, because it can produce continuous gradients. This article explores the use of B-splines to increase the accuracy of the multilevel summation method (for nonperiodic boundaries) without incurring additional computation other than a preprocessing step (whose cost also scales linearly). To obtain accurate results efficiently involves technical difficulties, which are overcome by a novel preprocessing algorithm. Numerical experiments demonstrate that the resulting method offers substantial improvements in accuracy and that its performance is competitive with an implementation of the fast multipole method in general and markedly better for Hamiltonian formulations of molecular dynamics. The improvement is great enough to establish multilevel summation as a serious contender for calculating pairwise interactions in molecular dynamics simulations. In particular, the method appears to be uniquely capable for molecular dynamics in two situations, nonperiodic boundary conditions and massively parallel computation, where the fast Fourier transform employed in the particle-mesh Ewald method falls short.

Born-Oppenheimer Ab Initio QM/MM Molecular Dynamics Simulations of Enzyme Reactions.

There are two key requirements for reliably simulating enzyme reactions: one is a reasonably accurate potential energy surface to describe the bond-forming/breaking process as well as to adequately model the heterogeneous enzyme environment; the other is to perform extensive sampling since an enzyme system consists of at least thousands of atoms and its energy landscape is very complex. One attractive approach to meet both daunting tasks is Born-Oppenheimer ab initio QM/MM molecular dynamics (aiQM/MM-MD) simulation with umbrella sampling. In this chapter, we describe our recently developed pseudobond Q-Chem-Amber interface, which employs a combined electrostatic-mechanical embedding scheme with periodic boundary condition and the particle mesh Ewald method for long-range electrostatics interactions. In our implementation, Q-Chem and the sander module of Amber are combined at the source code level without using system calls, and all necessary data communications between QM and MM calculations are achieved via computer memory. We demonstrate the applicability of this pseudobond Q-Chem-Amber interface by presenting two examples, one reaction in aqueous solution and one enzyme reaction. Finally, we describe our established aiQM/MM-MD enzyme simulation protocol, which has been successfully applied to study more than a dozen enzymes.

Spectral Ewald Acceleration of Stokesian Dynamics for polydisperse suspensions

In this work we develop the Spectral Ewald Accelerated Stokesian Dynamics (SEASD), a novel computational method for dynamic simulations of polydisperse colloidal suspensions with full hydrodynamic interactions. SEASD is based on the framework of Stokesian Dynamics (SD) with extension to compressible solvents, and uses the Spectral Ewald (SE) method [Lindbo and Tornberg (2010) [29]] for the wave-space mobility computation. To meet the performance requirement of dynamic simulations, we use Graphic Processing Units (GPU) to evaluate the suspension mobility, and achieve an order of magnitude speedup compared to a CPU implementation. For further speedup, we develop a novel far-field block-diagonal preconditioner to reduce the far-field evaluations in the iterative solver, and SEASD-nf, a polydisperse extension of the mean-field Brownian approximation of Banchio and Brady (2003) [39]. We extensively discuss implementation and parameter selection strategies in SEASD, and demonstrate the spectral accuracy in the mobility evaluation and the overall O(Nlog⁡N) computation scaling. We present three computational examples to further validate SEASD and SEASD-nf in monodisperse and bidisperse suspensions: the short-time transport properties, the equilibrium osmotic pressure and viscoelastic moduli, and the steady shear Brownian rheology. Our validation results show that the agreement between SEASD and SEASD-nf is satisfactory over a wide range of parameters, and also provide significant insight into the dynamics of polydisperse colloidal suspensions.