A generalized Ewald decomposition for screened Coulomb interactions

The interactions employed in molecular dynamics (MD) simulations are typically considered to be either short- or long-range. In d dimensions, interactions of the form 1/rⁿ are considered short range if $n>d$⁠. Short-range interactions are typically treated with a fast neighbor search within a cutoff radius r_c; the interaction is truncated thereafter, and a tail correction can be added.^1,2 Conversely, if $n≤d$⁠, then the interactions are considered long range; such interactions are considerably more difficult to treat. In periodic systems, most methods for handling long-range interactions are based on the Ewald sum,^3,4 which decomposes the Coulomb potential into two sums that are rapidly convergent. This decomposition can be implemented through a wide variety of algorithms, of which a naive implementation yields $O(N2)$ scaling; a simple optimization reduces this scaling to $O(N3/2)$⁠.⁵ Faster, grid-based methods have been developed that allow the fast Fourier transform (FFT) to be used,⁶ resulting in $O(NlogN)$ scaling; among these methods are the particle-mesh Ewald (PME),⁷ smoothed particle-mesh Ewald (SPME),⁸ and particle-particle particle-mesh (PPPM) methods.^6,9–12 Alternatively, non-Ewald methods¹³ have also been developed, such as the Wolf summation method¹⁴ and fast multipole methods,¹⁵ and hybrid methods^16,17 have been developed, as well.

However, this categorization of interactions as either short range or long range reflects the mathematical issue of force-sum convergence. In fact, many interactions can be considered to be medium range when their convergence is guaranteed yet either that there are a prohibitive number of particles in the neighbor radius or there are unacceptable truncation errors. For example, power-law potentials of the form $∼1/r2$ (charge-dipole), $∼1/r3$ (dipole-dipole), and 1/r⁶ (dispersion) fall into this category,^18,19 as well as Sutton-Chen²⁰ potentials relevant for metallic glasses.²¹ Isele-Holder et al.²² have extended the PPPM algorithm to interactions other than 1/r in their work on 1/r⁶ dispersion interactions. Similarly, prior work on screened Coulomb interactions includes extending the Ewald summation to Yukawa (exponentially screened) interactions.^23,24

Here, we examine a particular class of medium-range interactions that form with linearly screened Coulomb potentials, where the screening is specified through an arbitrary dielectric response function (DRF) $ϵ(𝐤)$⁠. The Ewald sum is formulated for this screened interaction and implemented with both a direct force approach (linked-cell list (LCL)²⁶ alone) and the PPPM method. An error analysis is carried out to examine the crossover between the direct force calculation and the PPPM implementation as a function of the range of the interaction and the particle number.

This paper is organized as follows. In Sec. II, we describe a class of interactions for screened charged systems obtained through the Fourier representation of the dielectric properties of background species. Several choices for the screening function are given in terms of the dielectric response function $ϵ(𝐤)$⁠, and their ranges are discussed. In Sec. III, the Ewald decomposition is developed for charged particle systems interacting through an effective interaction involving $ϵ(𝐤)$⁠.

Sec. IV includes error estimates and timing studies for the Yukawa dielectric response function.

Finally, we apply our results to compute the dynamic structure factor for a screened Coulomb system and discuss the implications of the choice of $ϵ(𝐤)$ and the computational costs that result from the temperature and density dependence of the effective interaction.

Consider N point particles with charges {Q_i} (e.g., ions), which can be partial, in the presence of a background charge density −en_b (e.g., electrons). The electrostatic forces in such a system can be calculated from the Poisson equation

where $Φ(𝐫)$ is the total electrostatic potential, $δ(𝐫)$ is the delta function, and {r_i(t)} are the particle coordinates at time t. Once this equation is solved, the corresponding force on the $ith$ particle is calculated to be $𝐅i=−Qi∇Φ(𝐫=𝐫i)$⁠. Expression (1) can be equivalently represented in Fourier space as

By assuming that the background density evolves instantaneously with the point particles, which is a valid assumption for a very light species, such as electrons, and approximating the structure of the background density as a linear response to the potential fluctuations induced by the presence of the point particles, we can express the background density in terms of a response function $χ(𝐤)$ as

We can write the electrostatic potential as

where the reciprocal of the DRF can be expressed in terms of the response function as

Finally, the force on the $ith$ particle can be represented in terms of the pair interactions as

where the effective interaction u_ij(r) is given in Fourier space as

When the background does not respond, $χ(𝐤)=0$ in (3). Such a completely uniform background corresponds to a dielectric response of unity in (5), and inverting (7) recovers the usual Coulomb interaction u_ij(r) = Q_iQ_j/r. Clearly, in the limit $ϵ→1$⁠, Ewald-class methods are needed. To continue further, we must make a specific choice of the DRF; we begin with the most widely used form

where $λs$ is the screening length associated with the background fluctuations. The associated pair interaction is given in real space by

This interaction decays exponentially and is therefore “short” range. However, as $λs→0$⁠, the number of particles in a neighbor-sphere diverges. It is therefore convenient to characterize the strength of the screening in terms of the dimensionless quantity

where $ai=(3/4πni)1/3$ is the ion-sphere radius, which is a measure of the interparticle distance. The $κ=0$ case is the pure Coulomb (one-component plasma) limit; thus, despite the exponential decay for finite κ, (9) has a long-range interaction as one of its limits. This suggests that there may be a finite range of κ for which the interaction (9) is medium range in the sense that there is no clear, optimal choice of algorithm.

To provide a concrete example, we can examine this choice in the context of the well-known Thomas-Fermi (TF) model,^27,28 for which

where T is the background temperature in energy units, E_F is the Fermi energy of the electron background in energy units, e is the elementary charge, and n_b is the background number density. Let us suppose that, for example, an acceptable error is given by a choice of cutoff radius $rc=5λs$⁠, which yields the factor $e−5∼0.674⋅10−2$ in (9). Allowing physical conditions (T and n_b) to vary, we can examine the number of particles in the neighbor cell. This is shown in Fig. 1, where contours show the logarithm of the number of particles in such a sphere as a function of temperature and density. Condensed systems (e.g., ionic liquids, liquid metals, and electrolytic solutions) tend to lie near the lower right region of the graph. Laser-produced plasmas originate as cold solids, in that same region, but can be rapidly heated to the very high temperatures of the upper right region. Roughly speaking, the lower right half of the diagram corresponds to $∼103$ particles in the neighbor cell for this choice of error. In Fig. 2, we show the $Γ−κ$ phase diagram. The coupling parameter Γ is given by

where Z is the degree of ionization and T_i is the ionic temperature. This diagram reflects the regimes of different charged particle systems of interest such as liquid metals,²⁹ ultracold neutral plasmas,³⁰ dusty plasmas,³¹ warm dense matter,³² and conditions at the National Ignition Facility (NIF).³³ We also show the regions corresponding to different liquid-like properties based on the properties of their velocity autocorrelation function.²⁷ Note that these disparate systems cover a wide range of $κ$ values.

Of course, the long-wavelength form of (8) is not unique, and interactions other than (9) are possible.³⁴ Moreover, the method is applicable to a system of partial charges as well. For this reason, we will retain the generality of an arbitrary $ϵ(𝐤)$ throughout the Ewald development in Sec. III. In fact, medium range interactions of different nature occur in biological and chemical systems where the linear screening approach is not applicable.^35–37 

To proceed with the development of the Ewald composition of (4), we note that the central quantity is the potential-like quantity

As is common, we implement the Ewald decomposition by adding and subtracting a screening-charge distribution of the form $ζρ(𝐫)$⁠, where $ζ$ is a constant to be determined, and $ρ(𝐫)$ is chosen to have the Gaussian form

where $α$ is the Ewald splitting parameter. Fourier-space representation of the decomposition takes the form

where the real-space contribution

can be highly short range (because there are effectively two sources of cutoff). Analogously,

is short range in Fourier-space and hence can also be computed efficiently; this is the subject of Sec. IV. In the absence of screening, $ϵ(k)→1$⁠, and we recover the Ewald decomposition for the standard Coulomb interaction.⁴ 

The real-space representation of $ϕR(𝐫)$ is given by the inverse Fourier transform

where an isotropic dielectric function $ϵ(k)$ is assumed in the second line. We evaluate this first with the long-wavelength (⁠$k→0$⁠) form of $ϵ(k)$ given by

where k_s is the magnitude of the inverse screening vector given by $ks=1/λs$⁠, where $λs$ is the long-wavelength screening length associated with properties of the background (i.e., classical, quantum, and multispecies). This form yields what is generically referred to as the Yukawa interaction. In what follows, we will be using $ks−1$ in place of $λs$⁠. Substituting $ϵY−1(k)$ in (20) yields

with $ζeks2/4α2=1$ (see Appendix  A). Alternative derivations of $ϕRY(r)$ have been given previously.^23–25 

Gradient corrections to the Yukawa result can be captured by the form

where leading-order corrections to the Yukawa form appear as k⁴ terms. An example of this type of DRF is the exact gradient-corrected screening (EGS) potential,³⁴ in which k⁴ terms arise from either electronic correlations or Heisenberg uncertainty (see Appendix  A). Interestingly, the form of the potential is qualitatively different for $ν<1$ and $ν>1$⁠, with the latter corresponding to an oscillatory decay.³⁴ Hence, the Ewald decomposition for the two cases will be different; in this paper, we will restrict our discussion to the $ν<1$ case. Substituting $ϵEGS−1(k)$ in (20) yields

where

Because $k−<k+$ (for $ν<1$⁠), the optimal choice is $ζ=e−k−2/4α2$⁠. We have given the detailed derivation in Appendix  A.

The real-space part of the Ewald decomposition of the screened Coulomb interaction between two point charges for an arbitrary DRF is, therefore, given by

This completes the formulation of the real-space contribution, which will be common to all Ewald-based methods. In Sec. IV, we first turn to the Fourier-space contribution, which we will treat with the PPPM algorithm for the particle-mesh contribution. Next, we provide an error analysis and a time-per-time-step comparison of the force calculations using PPPM and LCL algorithms for the Yukawa interaction. For the remainder of this paper, when we discuss LCL for a Yukawa interaction, this should be understood to refer to a force calculation using the LCL algorithm for a plain Yukawa interaction that is not subject to Ewald decomposition.

In this section, we begin with a brief explanation of the PPPM algorithm, followed by an error analysis for the computed force for the Yukawa DRF, with errors due to the truncation and discretization that are inherent to the algorithm. Then, we compare the execution times of the PPPM and LCL algorithms for the Yukawa potential to distinguish between the two algorithms in terms of computational efficiency for a certain level of error in the force calculation. The results of this section and Sec. V were obtained using the codes we implemented.³⁸ 

The particle-particle (PP) part of the PPPM algorithm refers to the real-space part of the Ewald decomposition, which can be highly short range; therefore, interactions beyond a certain cutoff radius can be neglected. If the number of particles within a cutoff sphere Nc is small compared to the total number of particles N in a simulation box, then the total number of interactions $O(NNc)$ tends to $O(N)$ as Nc becomes much smaller than N (⁠$Nc≪N$⁠). Such an $O(N)$ scaling is achieved by the LCL algorithm,⁶ in which the simulation box is divided into a number of cells with dimensions proportional to the cutoff radius. The particle-mesh (PM) portion of the PPPM algorithm refers to the Fourier-space part of the interaction that is computed on a grid using one of the two algorithms in Ref. 9; these algorithms differ in the manner in which forces are computed on a grid. Both versions of the PPPM algorithm employ the highly efficient FFT, which scales as $O(Mlog⁡M)$⁠, where M is the number of mesh points in the Fourier-space grid.

Errors are associated with computing the energy and forces in real-space and Fourier-space. Because we are interested in MD simulations in which the trajectories of a system of particles evolve primarily because of the forces acting on the particles, we restrict our error analysis to the errors in the forces. We chose to quantify the error in a force using the root-mean-square (RMS) error.^9,10 The RMS errors in the real-space and Fourier-space components of the forces, which are acting on particles in a simulation box, have a common form given by

where $P=R,F$ refers to real-space or Fourier-space, $𝐅i,Papprox.$ is the approximate force acting on particle i, $𝐅i,Pexact$ is the actual force (or numerically exact force) acting on particle i, and $ξ$ is used to distinguish between computed RMS errors and analytical estimates of RMS errors. Moreover, because the computational costs associated with real-space and Fourier-space force calculations are determined by the levels of error that are allowed, it is essential to understand how the errors vary with respect to the parameters. For this purpose, we have derived analytical estimates for the errors in the real-space and Fourier-space components of the forces for a system of randomly distributed particles. The error estimates were verified by comparing with the RMS errors computed for a system of particles placed at random positions within a simulation box. For the remainder of this paper, we use a_i as the unit of distance and the elementary charge e as the unit of charge; thus, the RMS force error is in units of $e2/ai2$⁠.

The error in the real-space component of the force is due to the truncation of particle interactions beyond a certain cutoff radius. For a system of N charged point particles interacting pair-wise with an arbitrary short-range interaction expressed as $us(𝐫i−𝐫j)=QiQjϕs(𝐫i−𝐫j)$⁠, the estimate of the RMS force error due to a cutoff radius r_c can be derived by adapting the error analysis in Ref. 39, which considers the real-space portion of the Ewald decomposition for the Coulomb potential. With the assumption that beyond r_c particles are randomly distributed, the estimate of the RMS force error due to u_s(r) is given by

where $𝐟s(𝐫)=−∇𝐫ϕs(𝐫)$⁠, $Q=∑i=1NQi2$⁠, and V is the volume enclosing the particles.

Substituting $𝐟RY(𝐫)=−∇𝐫ϕRY(𝐫)$ in (29) and making several approximations, we get the following estimate of the RMS force error for the real-space portion of the Ewald decomposition for the Yukawa potential:

(see Appendix  B for details). Note that as $κ→0$ (Coulomb limit), we recover the corresponding error estimate for the Coulomb case.³⁹ Figs. 3(a)–3(c) show very good agreement between the error estimates and the computed RMS error in the real-space force component for a system of $5×103$ unit charges with $κ$ = 0.1, 1, and 2 and varying values of $α$⁠. The choice of these $κ$ values is consistent with the $κ$ range of interest as shown in Fig. 2. The particles were randomly distributed in a cube of edge length L given by $L=(4πN/3)1/3$ in units of a_i.

In the Fourier-space part of the PPPM algorithm, the force calculation for a system of charged point particles with arbitrary spatial distribution is mapped to a real-space grid by assigning particle positions to the real-space grid points using B-splines of some order p. Calculations are then performed in Fourier-space starting with the FFT of the charge density. This is followed by computing the potential energy and forces on the Fourier-space grid following the steps described in Ref. 9. Forces on the Fourier-space grid are then transformed to real-space by an inverse FFT. Next, the forces are interpolated from the real-space grid to the actual particle positions using B-splines of the same order, in a manner similar to the earlier step of assigning particle positions to the real-space grid. FFT involves an inherent truncation in the maximum allowable length of Fourier-space vectors, which is one cause for error in the Fourier-space component of the forces. In addition, errors occur because of aliasing, which arises because calculations are being made on a grid, and because of the interpolation of particle positions to the real-space grid and of forces computed on the real-space grid to actual particle positions. Critically, the PPPM algorithm employs an optimal Green function⁶ that is derived by minimizing the error in the forces for a random distribution of particles. The optimal Green function can be interpreted as a weighted average of Fourier modes of the actual Green function.

For the version of the PPPM algorithm in which forces are computed on a grid,⁹ the optimal Green function corresponding to an arbitrary Green function is given by⁹ 

where n refers to the grid vector denoted by the triplet of grid indices (n_x, n_y, n_z), $G^𝐤𝐧$ is the actual Green function evaluated at the Fourier-space vector k_n corresponding to grid vector n, m refers to the triplet of grid indices (m_x, m_y, m_z) that contribute to aliasing, and $U^𝐤𝐧$ is the Fourier transform of the B-spline of order p. $U^𝐤𝐧$ is given by

where M_x, M_y, and M_z refer to the number of grid points along the respective directions in Fourier-space. The triplet of grid sizes (h_x, h_y, h_z) are then given by (L_x/M_x, L_y/M_y, L_z/M_z). The estimate of the optimal RMS error in the Fourier-space force is given by⁹ 

where

Figs. 3(d)–3(f) show very good agreement between the RMS error estimates for the Fourier-space force with order of B-splines ranging from 2 to 7 and the corresponding RMS error computed using the same system of $5×103$ unit charges that was employed for the real-space error analysis. Because the simulation volume is a cube, the numbers of grid points M along each direction in Fourier-space is equal.

Eq. (33) is not of a form that explicitly reveals the influence of the parameters h, α and p on the error incurred in the Fourier-space force calculation. Therefore, (33) requires further simplification, which is possible if the arbitrary Green function decays sufficiently rapidly in Fourier-space that the alias contributions can be neglected. Following the approach in Ref. 10, (33) is approximated to

where $AF$ is given by

where $β(p,m)$ denotes an integral given by

and $Cm(p)$ are coefficients listed in Table I of Ref. 10 (see Appendix  B for further detail about the simplification of (33)). For the Yukawa interaction, substituting $G^k=4π(k2+κ2)e−(k2+κ2)/4α2$ results in a modified form for $β(p,m)$ given by

where $ψ(p,m,κ/α)$ denotes an integral of the form

Finally, a simplified form of the RMS error estimate for the Fourier-space component of the force for the Yukawa interaction is given by

where

Fig. 4 shows a good agreement between $ΔFFY(approx.)$ and the full error estimate $ΔFFY$ for $κ$ = 0.1 and 1 for a system of 10³ unit charges. For $κ=1$ as $α$ decreases below unity, $ΔFFY(approx.)$ deviates from $ΔFFY$⁠. Further, we found that for $κ>1$⁠, $ΔFFY(approx.)$ is not a good approximation of $ΔFFY$ for the $α$ values of interest. Nevertheless, the simplified form of $ΔFFY(approx.)$ is useful because it reveals that the error approximately varies as the product of the grid size and the Ewald splitting parameter (⁠$hα$⁠) raised to the power given by the order of B-spline (p). Further, finding an optimal set of parameters requires considerable effort using (33), whereas using a simplified estimator can expedite the process of finding the optimal parameters. An alternative approach for estimating the errors is described in Ref. 54.

The total RMS error in the force computed using the PPPM algorithm is defined as

To compare the computational costs of the PPPM and LCL algorithms, the execution times for the two algorithms must be compared for the same RMS force errors. To this end, an estimate of the RMS force error for Yukawa interaction is needed. Substituting $𝐟Y(𝐫)=−∇𝐫ϕY(𝐫)$ in (29), where $ϕY(𝐫)=e−κr/r$⁠, results in the following estimate of the RMS force error for Yukawa interaction:

This equation, in turn, implies that

(see Appendix  B for further detail). As shown in Fig. 5, the RMS error estimate in (43) has been verified by comparison with the RMS error computed using 10⁵ unit charges placed at random positions in a box.

Based on the error analysis for the Yukawa interaction, Fig. 6 shows that for $κ=0.1$ and an RMS force error of 10⁻⁸ (⁠$e2/ai2$⁠), the minimum image convention $rc≤L/2$ is not satisfied below a system size of 10⁷ unit charges. Therefore, for a $κ$ value and particle numbers for which LCL would be less efficient, PPPM serves as an excellent alternative. As $κ$ increasesfor approximately the same error of 10⁻⁸ (⁠$e2/ai2$⁠), the threshold above which the particle number satisfies the minimum image convention decreases; for example, for $κ=0.5$⁠, the minimum image convention is satisfied beyond 10⁵ unit charges, as shown in Fig. 6. Importantly, for some conditions where the LCL algorithm is applicable, we find the PPPM algorithm to be computationally advantageous. This observation is based on the comparison of the time per time step of PPPM and LCL algorithms for a fixed RMS force error.The timing results for a range of $κ$ and N (with unit charges) are provided in Table I. For example, with 10⁶ unit charges and an RMS force error of approximately 10⁻⁵ (⁠$e2/ai2$⁠), we find that the PPPM algorithm was approximately 25 and 130 times faster than the LCL algorithm for $κ=0.5$ and 0.3, respectively (as shown in Table I). As $κ$ approaches unity, we find the LCL and PPPM algorithms to be comparable in performance, as summarized by the entries in Table I for $κ=1$⁠. As $κ$ increases beyond unity, we find that the PPPM algorithm loses its computational advantage over the LCL algorithm, as can be seen in the entries in Table I for $κ=2$⁠. We should point out, however, that execution times are architecture dependent and our conclusions are based on that fact. All calculations used a 3.5 GHz Intel Xeon(R) processor and 32 GB RAM.

In this section, we apply our results to the computation of the dynamic structure factor $S(k,ω)$⁠. Our interest in $S(k,ω)$ falls into two categories. Theoretically, $S(k,ω)$ contains all of the many-body processes associated with thermal density fluctuations, including collective modes and transport.^29,40–42 All of these processes are obviously impacted by the choice of $ϵ(k)$⁠, and we wish to use MD for model validation. The ionic $S(k,ω)$ is also connected directly with neutron scattering⁴³ and, indirectly, to X-ray Thomson Scattering (XRTS).⁴⁴ In the context of warm dense matter physics, for example, recent XRTS studies have developed a modified Navier-Stokes equation (MNSE) to model the ionic $S(k,ω)$⁠.²⁹ 

Within the framework of the Yukawa model, we are interested in density and temperature conditions for which MD is required to validate approximate theoretical models and where LCL would be computationally less efficient than PPPM for force calculations. Consider the case $Γ=2$ and $κ=0.1$⁠, for 10⁴ unit charges the PPPM algorithm can have a force error of approximately 10⁻⁶ (⁠$e2/ai2$⁠) for an optimal choice of $α$⁠, r_c, h, and p. For the same error and particle number, LCL would be less efficient because the minimum-image convention is not satisfied. This case corresponds to an effective coupling of $Γeff=Γe−κ≃0.74$⁠, which is of an order that is relevant to systems of recent interest, such as warm dense matter.⁴⁴ 

The dynamic structure factor for an N-particle system is defined as

Here, the brackets $⟨⋯⟩$ refer to the stationary time average

where T is the length of the run. We note that for a system in equilibrium the stationary time average is equivalent to the ensemble average as $T→∞$⁠. The Fourier-transformed densities n(k,t) are obtained from the definition of the density,

to yield

Substituting (46) into (45), we obtain

which, through the change of variables $y=t+τ$⁠, can be written as

We use the usual definition of the temporal Fourier transform

For numerical purposes, the infinite domain must be approximated as finite, and as we are calculating stationary processes, we can take this domain to be [0, T] without loss of generality for sufficiently large T. This yields the approximate transform

Therefore, the final form for $S(𝐤,ω)$ is given by

as $n(−𝐤,−ω)$ = $n*(𝐤,ω)$⁠. The Fourier-transformed density $n(𝐤,ω)$ is given by

Particle positions from MD are available only at certain times in [0, T]; hence, $n(𝐤,ω)$ is computed using

where N_t is the number of steps after the equilibration phase of MD, and $Δt$ is the step size for time integration. Particle positions are stored at intervals corresponding to this step size. Long runs are preferred because the noise in $S(𝐤,ω)$ decreases as T increases. The temporal Fourier transform is computed by FFT. $S(𝐤,ω)$ computed from MD data tends to be noisy, and the extent of the noise increases as k increases. To determine the location of the peak of $S(𝐤,ω)$ and its width, the noisy curves are usually smoothed by averaging over data from different MD runs that differ in their initial conditions. In Sec. V B, we present a potentially more efficient smoothing approach that employs kernel smoothing.^49,50

Kernel smoothing refers to weighted averaging of noisy observations to obtain a smooth function.^49,50 To smooth our data using kernel smoothing, we use a kernel regression technique⁵⁰ defined by

where y(x_i) is the observation of y at x_i, B is the number of observations, and $Kσ(x∼,xi)$ is the smoothing kernel. This yields a weighted average to give a smooth function $y∼(x∼)$ evaluated at $x∼$⁠. We use the Gaussian kernel

where $σ$ is the kernel width. For $Γ=50$ and $κ=1$⁠, Figs. 7 and 8 show the smooth dynamic structure factor curves (red curves) obtained by kernel smoothing of noisy dynamic structure factor curves obtained from MD simulations (grey curves) with σ = 0.01 and σ = 0.005, respectively. MD simulations of 10³ unit charges were performed using the PPPM algorithm, with an RMS force error of approximately 10⁻⁴ (⁠$e2/ai2$⁠). The parameters used were $α$=1, r_c=3, M=64, and p = 2.

We have defined a simple way of quantifying the smoothing error as follows:

where $y∼2(x∼)$ is defined as

More rigorous ways of quantifying the error in smoothing can be found in Refs. 49 and 50. Figs. 7 and 8 show the error in smoothing (blue band) about the smooth $S(q,Ω)$ curves (red curves). Because the kernel width influences the extent of smoothing and the corresponding error in the smoothing, the width must be chosen carefully such that the smoothed data are consistent with the physically expected result.

In this section, we compare dispersion relations obtained using PPPM-MD with those obtained using different theoretical models. For convenience, k and ω are expressed as dimensionless quantities q = ka_i and $Ω=ω/ωp$⁠, where $ωp=(4πniQi2/M)1/2$ is the ion plasma frequency for ionic mass M and ionic number density n_i.

One simple dispersion relation obtained from a static structure factor S(q)²⁹ is given by

In the random-phase approximation (RPA) for a Yukawa system, S(q) takes the approximate form

More accurate theoretical forms of S(q) are discussed in Refs. 29 and 46. Fig. 9(a) shows the comparison of a dispersion relation obtained using PPPM-MD and dispersion relations obtained using different theoretical approximations. The purpose of this comparison is to show that the MD result agrees well with the theoretical dispersion relations in the low-q limit, where the theory is considered to work well.²⁹ As q increases, the MD curve begins to deviate from the theoretical curves; this occurs because of the lack of required many-body effects in those theoretical models. A very accurate S(q) can be obtained from hypernetted chain (HNC) equations with a bridge function⁴⁸ (denoted by S^HNC(q)). Although particle correlations are captured in $Ω(q)$ by S^HNC(q) (solid blue curve), this model does not incorporate viscoelastic effects,⁴⁶ and hence it underestimates $Ω(q)$ compared to MD. Fig. 9(b) shows the full-width-half-maximum (FWHM) of a smooth $S(q,Ω)$ constructed by kernel smoothing of $S(q,Ω)$ obtained from MD data using a Gaussian kernel of widths $σ=$0.01 (blue dots) and $σ=$ 0.005 (green dots). A smaller width of $σ=$0.005 results in a curve that approaches zero (as expected from theory) more accurately than does a curve constructed using a width of $σ=$ 0.01.