Large Number Of Grid Points Research Articles

Superpositions and higher order Gaussian beams

High frequency solutions to partial differential equations (PDEs) are notoriously difficult to simulate numerically due to the large number of grid points required to resolve the wave oscillations. In applications, one often must rely on approximate solution methods to describe the wave field in this regime. Gaussian beams are asymptotically valid high frequency solutions concen- trated on a single curve through the domain. We show that one can form integral superpositions of such Gaussian beams to generate more general high frequency solutions to PDEs. As a particular example, we look at high frequency solutions to the constant coefficient wave equation and construct Gaussian beam solutions with Taylor expansions of several orders. Since this PDE can be solved via a Fourier transform, we use the Fourier transform solution to gauge the error of the constructed Gaussian beam superposition solutions. Furthermore, we look at an example for which the solution exhibits a cusp caustic and investigate the order of magnitude of the wave amplitude as a function of frequency at the tip of the cusp. We show that the observed behavior is in agreement with the predictions of Maslov theory.

Dynamic dipole polarizabilities of the ground and excited states of confined hydrogen atom computed by means of a mapped Fourier grid method

AbstractAccurate theoretical estimates of static and dynamic dipole polarizabilities are reported for the ground and excited states of hydrogen atom confined at the center of a spherical “box” with impenetrable walls using a novel theoretical algorithm that combines the variational–perturbation approach with an appropriately adapted mapped Fourier grid method and uniformly maintains the numerical accuracy. A variation of computed polarizabilities is observed as a function of the number of grid points. However, rapid convergence to their correct values, even far into the anomalous dispersion region, is achieved by an extrapolation procedure to the limit of an infinitely large number of grid points using a small number of the lowest‐order Padé approximants. It is shown that dipole polarizabilities strongly depend upon the electronic state and the radius of confinement. In particular, the static polarizability of 2s state changes sign under strong confinement. © 2007 Wiley Periodicals, Inc. Int J Quantum Chem, 2008

Dynamically adapted grids for interacting discontinuous solutions

Further development of the dynamic adaptation method for gas dynamics problems that describe multiple interactions of shock waves, rarefaction waves, and contact discontinuities is consid- ered. Using the Woodward-Colella problem and a nonuniformly accelerating piston as examples, the efficiency of the proposed method is demonstrated for the gas dynamics problems with shock wave and contact discontinuity tracking. The grid points are distributed under the control of the diffusion approx- imation. The choice of the diffusion coefficient for obtaining both quasi-uniform and strongly nonuni- form grids for each subdomain of the solution is validated. The interaction between discontinuities is resolved using the Riemann problem for an arbitrary discontinuity. Application of the dynamic adapta- tion method to the Woodward-Colella problem made it possible to obtain a solution on a grid consisting of 420 cells that is almost identical to the solution obtained using the WENO5m method on a grid con- sisting of 12 800 cells. In the problem for a nonuniformly accelerating piston, a proper choice of the dif- fusion coefficient in the transformation functions makes it possible to generate strongly nonuniform grids, which are used to simulate the interaction of a series of shock waves using shock wave and contact discontinuity tracking. The dynamic adaptation methods (13-21) are less complicated and make it possible to generate grids with a constant (19) or variable (18) number of grid points. In these methods, the control of the distribution of the grid points is achieved by using information about the dynamic behavior of the solution to be found. This makes it possible to concentrate a large number of grid points in the regions where the solution under- goes sharp changes (13). A close relationship between the dynamic behavior of the solution and the location of the grid points requires that the coordinates of the grid points be recalculated at each time layer. This fact places high requirements upon the consistency of the solution dynamics with the grid points motion. For this reason, the algorithms that do not use any fitting parameters have certain advantages (15). These fea- tures are especially pronounced in the generation of adaptive grids for unsteady gas dynamics problems (see (16-22)), which describe rapidly changing processes. The solution of hyperbolic differential equations,

A Boltzmann-based finite volume algorithm for surface water flows on cells of arbitrary shapes

An explicit two–dimensional conservative finite volume model for shallow water equations is formulated and tested. The algorithm for the mass and momentum fluxes at the control surface of the finite volume is obtained from the solution of the Bhatnagar–Gross–Krook (BGK) Boltzmann equation. Unlike classical methods, BGK schemes do not require an ad–hoc splitting of advection and diffusion. The BGK scheme is second order in both time and space. The formulation of the BGK algorithm is performed for a cell of arbitrary irregular shape, but the test cases are conducted using a structured grid of quadrilateral cells. Two approximate Riemann solvers, the HLLC scheme and the two–stage Hancock–HLLC scheme, where HLL stands for Harten, Lax and van Leer and C stands for contact discontinuity, are also considered. The second–order accuracy of HLL and Hancock–HLLC schemes is obtained by MUSCL approach, where MUSCL is the acronym for Monotone Upstream–centered Schemes for Conservation Laws. The data reconstruction for all three schemes is carried out by theVan Leer limiter. The test cases involve strong shocks and expansion waves. The accuracy of the schemes are measured using an absolute error norm and a waviness error norm. The HLLC scheme is highly oscillatory for Courant number larger than 0.5, while the BGK and the Hancock–HLLC schemes are applicable for Courant numbers as high as 1.0. For a fixed value of the central processing unit (CPU) time, the absolute error of the Hancock–HLLC is slightly smaller than that of the BGK while the waviness error of the BGK is quite close to that of Hancock–HLLC. This is because (i) the Hancock–HLLC is a two–step method while the BGK is a single–step method (i.e., the Hancock–HLLC requires storage of intermediate variables, but the BGK does not), and (ii) the Hancock–HLLC schemes requires larger number of grid points than the BGK scheme for the same level of accuracy. For example, to achieve an absolute error of 0.01, the BGK requires about 600 grid points while the Hancock–HLLC requires about 800 grid points. Both the BGK and Hancock–HLLC schemes have similar convergence properties. Unlike exact or approximate Riemann solvers, BGK fluxes accounts for both waves and diffusion. The ability of the BGK scheme to model diffusion is illustrated using a viscous flow problem. Excellent agreement between the analytical and computed viscous flow solution is found. Although the BGK and Hancock–HLLC schemes perform similarly for hyperbolic problems, BGK schemes have the added advantage of being able to solve hyperbolic–parabolic problems without the need for an ad–hoc operator splitting. This is important given that the artificial splitting of advection and diffusion is known to cause artificial widening in shear layers and introduces artificial transient in regions with sharp gradients. Such problems arise when the splitting operation fails to faithfully represent the correct coupling between the physics of advection and the physics of waves.

Lobatto Point Quadrature for Thermal Lubrication Problems Involving Compressible Lubricants. EHL Applications

Refined solutions of thermal lubrication problems generally require fine mesh and many iteration steps. To resolve these difficulties, Elrod and Brewe proposed an efficient algorithm based on the use of a Lobatto point version of the Gauss quadrature, which is typically twice as fast as the other quadrature methods. The original Lobatto algorithm was only applied for incompressible hydrodynamic lubrication. This paper presents a Lobatto point quadrature algorithm which is applicable for thermal elastohydrodynamic lubrication (EHL) problems where both density and viscosity of the lubricant are taken to be temperature and pressure dependent and the transverse velocity term in the energy equation is obtained from the continuity equation. Within this approach, the unknown temperature across the film is written in a series of Legendre polynomials. Regardless of the order of the series expansions, the thermal Reynolds equation can explicitly contain only the information from the first three Legendre polynomials, i.e., data from up to a second-order polynomial. Use of the Lobatto point calculation method has resulted in accuracy without the use of a larger number of grid points.

A Fast Parallelized Computational Approach Based on Sparse LU Factorization for Predictions of Spatial and Time-Dependent Currents and Voltages in Full-Body Biomodels

Realistic and accurate numerical simulations of electrostimulation of tissues and full-body biomodels have been developed and implemented. Typically, whole-body systems are very complex and consist of a multitude of tissues, organs, and subcomponents with diverse properties. From an electrical standpoint, these can be characterized in terms of separate conductivities and permittivities. Accuracy demands good spatial resolution; thus, the overall tissue/animal models need to be discretized into a fine-grained mesh. This can lead to a large number of grid points (especially for a three-dimensional entity) and can place prohibitive requirements of memory storage and execution times on computing machines. Here, the authors include a simple yet fast and efficient numerical implementation. It is based on LU decomposition for execution on a cluster of computers running in parallel with distributed storage of the data in a sparse format. In this paper, the details of electrical tissue representation, the fast algorithm, the relevant biomodels, and specific applications to whole-animal studies of electrostimulation are discussed

Ab initioinvestigation of noncontact atomic force microscopy tip-surface instability in InAs(110) surface

In the present work we report ab initio pseudopotential calculations based on density functional theory to investigate the noncontact atomic force microscopy (NC-AFM) image contrast on the InAs(110) $(1\ifmmode\times\else\texttimes\fi{}1)$ surface. The foremost tip structure is modeled by a $\mathrm{Si}{\mathrm{H}}_{3}$ tip. The effect of the tip-induced surface relaxations on the calculated forces was investigated for the tip above As and In atoms. The force curves corresponding to these vertical scans show an hysteretic behavior and this effect causes an energy dissipation of 0.3 (tip on top of As) and $1.8\phantom{\rule{0.3em}{0ex}}\mathrm{eV}$ (tip on top of In), respectively. The presence of this hysteresis suggests that stable NC-AFM images can be obtained for tip-sample distances before this instability. In this stable regime the force curves obtained for perturbed (due to tip-sample interaction) and unperturbed InAs(110) surface exhibit the same qualitative behavior. From the calculated forces for the unperturbed InAs(110) surface on a large number of grid points in real space we obtained maps of constant frequency shifts. The influence of long-range van der Waals forces on the simulated AFM images due to the macroscopic part of the tip was taken into account by an empirical model. The overall structure and the corrugation of the simulated NC-AFM images are in good agreement with the experimental results and allow us to explain the experimentally observed features of the image contrast mechanism on the basis of the calculated short-range chemical tip-sample interaction forces.

A New Efficient Algorithm for Solving an Incompressible Flow on Relatively Fine Mesh

The new coupled solver SOAR (SIMPLE Optimized and Automated Relaxation factors) proposed by the author has been reformulated on a collocated grid and compared with the SIMPLE method for computing performance as the grid is refined. The functional relation between calculational time T and total grid points N has been studied through 22 calculations for a rectangular tank and can be expressed as T ≈ k × Nm. Here, the value of the exponent m has been found to range from 1.2 to 1.5 for the SOAR solver and from 2.4 to 3.0 for the SIMPLE solver, and it is clear that the SOAR method has better computing performance than the SIMPLE method as the grid is refined. For large numbers of grid points, which are increasingly required for assessment of confidence in computational fluid dynamics, this translates into significant and predictable reductions in run times and faster project turnaround.

速度・応力型差分法での固体・流体境界の扱いについて

We must consider the seismo-acoustic problem at the irregular fluid-solid boundary in various acoustical and seismological subjects: e. g., the ocean-acoustics problem, the analysis of the data from ocean-bottom seismometers, the analysis of the scattered waves from the irregular core-mantle boundary, and the strong motion prediction at sites close to the ocean. In this paper, we study why and how the fluid-solid or the acoustic-seismic boundary conditions are satisfied in the velocity-stress staggered grid finite-difference method (FDM) which has been widely used for seismo-acoustics problems. A simple analysis of the 2D finite-difference equations together with the comparisons between finite-difference seismograms and the discrete wavenumber method (DWM) seismograms show that the conditions are satisfied provided that (1) the boundary is placed through the shear stress (τxz) grid points (i. e., the rigidity of the shear stress grid points on the boundary must be set to zero), (2) the 2nd-order centered equation with averaged density is applied to the normal velocity component of the grid points right on the boundary, and (3) inadequate higher order finite-difference equations (e. g. the standard 4th-order centered equation) is not applied to grid points neighboring the boundary so that the values with discontinuity at the boundary are not included in the finitedifference equations for grid points off from the boundary. The accuracy of the spatial finitedifference equation at the boundary is O(h) with h being the space increment despite its apparent form of “2nd-order” equation. Violations to these conditions deteriorate the results of computation, especially for the interface wave such as Stoneley wave. In FDM inclined interface is approximated by stairsteps made from the grid. We find that these stairsteps can generate scattered waves with amplitudes up to 20% of the Stoneley wave. We need about 60 grid points per wavelength of the Stoneley wave at the dominant frequency in order to suppress such scattered wave. In the case of basin-like fluid layer with free surface and with “shoreline” we need large number of grid points (about 120 grid points) per wavelength at the dominant frequency in order to model the secondary Rayleigh wave generated (scattered) at the basin edges.

Stochastic Strong Ground-Motion Simulation of the 7 September 1999 Athens (Greece) Earthquake

The stochastic method for finite faults is applied to simulate strong ground motion from the 7 September 1999, moment magnitude M 5.9 Athens earthquake. The method includes descritization of the fault plane into a certain number of subfaults, each of which is assigned an ω -2 spectrum. A slip-distribution model, derived from previous studies of this earthquake, is used to specifically account for the source effect. Contributions from all subfaults are then empirically attenuated to the observation sites, where they are summed to produce the synthetic acceleration time history. The method is first calibrated against its ability of reproducing the recordings at 19 strong-motion stations, at epicentral distances ranging from 16 to 61 km. The calibrated model is then applied to calculate synthetics at a large number of grid points covering the area around the fault plane. Simulated peak values are subsequently used to produce synthetic peak ground acceleration and spectral acceleration maps at hard rock. Both peak ground acceleration and spectral acceleration maps imply energy directivity toward the east, where most of the damage was concentrated. The directivity effect is most prominent at large periods (>2 sec) and in the period range 0.2 to 0.3 sec. Independent geotechnical studies showed considerable site effect at periods <0.5 sec within the meizoseismal area. This result, coupled with the results of the present study, imply that the damage distribution pattern of the 1999 Athens earthquake can be explained by the destructive combination of two factors: the source directivity and the site effect.

Comparison of discrete singular convolution and generalized differential quadrature for the vibration analysis of rectangular plates

This paper presents a comprehensive comparison study between the discrete singular convolution (DSC) and the well-known global method of generalized differential quadrature (GDQ) for vibration analysis so as to enhance the understanding of the DSC algorithm. The DSC method is implemented through Lagrange's delta sequence kernel (DSC-LK), which utilizes local Lagrange polynomials to calculate weighting coefficients, whereas, the GDQ requires global ones. Moreover, it is shown that the treatments of boundary conditions and the use of grid systems are different in the two methods. Comparison study is carried out on 21 rectangular plates of different combinations of simply supported (S), clamped (C) and transversely supported with nonuniform elastic rotational restraint (E) edges, and five rectangular plates of mixed supporting edges, some of which with a range of aspect ratios and rotational spring coefficients. All the results of the DSC-LK agree very well with both those in the literature and newly computed GDQ results. Furthermore, it is observed that the DSC-LK performs much better for plates vibrating at higher-order eigenfrequencies. Unlike the GDQ, the DSC-LK is numerically stable for problems which require a large number of grid points.

The State of the Art in Flow Visualisation: Feature Extraction and Tracking

Abstract Flow visualisation is an attractive topic in data visualisation, offering great challenges for research. Very large data sets must be processed, consisting of multivariate data at large numbers of grid points, often arranged in many time steps. Recently, the steadily increasing performance of computers again has become a driving force for new advances in flow visualisation, especially in techniques based on texturing, feature extraction, vector field clustering, and topology extraction. In this article we present the state of the art in feature‐based flow visualisation techniques. We will present numerous feature extraction techniques, categorised according to the type of feature. Next, feature tracking and event detection algorithms are discussed, for studying the evolution of features in time‐dependent data sets. Finally, various visualisation techniques are demonstrated. ACM CSS: I.3.8 Computer Graphics—applications

Time-domain boundary conditions for outdoor ground surfaces

Finite-difference time-domain techniques are promising for detailed dynamic simulations of sound propagation in complex atmospheric environments. Success of such simulations requires the development of new techniques to accurately handle the reflective and absorptive properties of a porous ground. One method of treating the ground boundary condition in the time domain [Salomons et al., Acta Acust. 88, 483–492 (2002)] is to use modified fluid dynamic equations, where the ground is considered as a porous medium described by its physical properties. However, this approach significantly increases computation time, as the domain must be extended into the ground and a large number of grid points are needed. Standard impedance models for the ground boundary condition are frequency-domain models, which generally are non-causal [Y. H. Berthelot, J. Acoust. Soc. Am. 109, 1736–1739 (2001)]. The development of a time-domain boundary condition from these models requires removing the singularity from the impedance equation when transforming from the frequency domain to the time domain. Alternatively, as the impedance boundary condition is a flux equation, a time-domain boundary condition can be derived from first principles, using the physical properties of the ground. We report on our development of a time-domain ground boundary condition.

On the influence of numerical schemes and subgrid–stress models on large eddy simulation of turbulent flow past a square cylinder

AbstractInfluence of finite difference schemes and subgrid‐stress models on the large eddy simulation calculation of turbulent flow around a bluff body of square cylinder at a laboratory Reynolds number, has been examined. It is found that the type and the order of accuracy of finite‐difference schemes and the subgrid‐stress model for satisfactory results are dependent on each other, and the grid resolution and the Reynolds number. Using computational grids manageable by workstation‐level computers, with which the near‐wall region of the separating boundary layer cannot be resolved, central‐difference schemes of realistic orders of accuracy, either fully conservative or non‐conservative, suffer stability problems. The upwind‐biased schemes of third order and the Smagorinsky eddy‐viscosity subgrid model can give reasonable results resolving much of the energy‐containing turbulent eddies in the boundary layers and in the wake and representing the subgrid stresses in most parts of the flow. Noticeable improvements can be obtained by either using higher order difference schemes, increasing the grid resolution and/or by implementing a dynamic subgrid stress model, but each at a cost of increased computational time. For further improvements, the very small‐scale eddies near the upstream corners and in the laminar sublayers need to be resolved but would require a substantially larger number of grid points that are out of the range of easily accessible computers. Copyright © 2002 John Wiley &amp; Sons, Ltd.

Multigrid Block Preconditioning for a Coupled System of Partial Differential Equations Modeling the Electrical Activity in the Heart

The electrical activity of the heart may be modeled with a system of partial differential equations (PDEs) known as the bidomain model. Computer simulations based on these equations may become a helpful tool to understand the relationship between changes in the electrical field and various heart diseases. Because of the rapid variations in the electrical field, sufficiently accurate simulations require a fine-scale discretization of the equations. For realistic geometries this leads to a large number of grid points and consequently large linear systems to be solved for each time step. In this paper, we present a fully coupled discretization of the bidomain model, leading to a block structured linear system. We take advantage of the block structure to construct an efficient preconditioner for the linear system, by combining multigrid with an operator splitting technique.

Parallel Algorithms with Local Fourier Basis

We present a nonoverlapping domain decomposition method with local Fourier basis applied to a model problem in liquid flames. The introduction of domain decomposition techniques in this paper is for numerical and parallel efficiency purposes when one requires a large number of grid points to catch complex structures. We obtain then a high-order accurate domain decomposition method that allows us to generalize our previous work on the use of local Fourier basis to solve combustion problems with nonperiodic boundary conditions (M. Garbey and D. Tromeur-Dervout, J. Comput. Phys.145, 316 (1998)). Local Fourier basis methodology fully uses the superposition principle to split the searched solution in a numerically computed part and an analytically computed part. Our present methodology generalizes the Israeli et al. (1993, J. Sci. Comput. 8, 135) method, which applies domain decomposition with local Fourier basis to the Helmholtz's problem. In the present work, several new difficulties occur. First, the problem is unsteady and nonlinear, which makes the periodic extension delicate to construct in terms of stability and accuracy. Second, we use a streamfunction biharmonic formulation of the incompressible Navier–Stokes equation in two space dimensions: The application of domain decomposition with local Fourier basis to a fourth-order operator is more difficult to achieve than for a second-order operator. A systematic investigation of the influence of the method's parameters on the accuracy is done. A detail parallel MIMD implementation is given. We give an a priori estimate that allows the relaxation of the communication between processors for the interface problem treatment. Results on nonquasi-planar complex frontal polymerization illustrate the capability of the method.

Scattering by hydraulic fractures: Finite‐difference modeling and laboratory data

Reservoir production can be stimulated by creating hydraulic fractures that effectively facilitate the inflow of hydrocarbons into a well. Considering the effectiveness and safety of the operation, it is desirable to monitor the size and location of the fracture. In this paper we investigate the possibilities of using seismic waves generated by active sources to characterize the fractures. First, we must understand the scattering of seismic waves by hydraulic fractures. For that purpose we use a finite‐difference modeling scheme. We argue that a mechanically open hydraulic fracture can be represented by a thin, fluid‐filled layer. The width or aperture of the fracture is often small compared to the seismic wavelength, which forces us to use a very fine grid spacing to define the fracture. Based on equidistant grids, this results in a large number of grid points and hence computationally expensive problems. We show that this problem can be overcome by allowing for a variation in grid spacing in the finite‐difference scheme to accommodate the large‐scale variation in such a model. Second, we show ultrasonic data of small‐scale hydraulic fracture experiments in the laboratory. At first sight it is difficult to unravel the interpretation of the various events measured. We use the results of the finite‐difference modeling to postulate various possible events that might be present in the data. By comparing the calculated arrival times of these events with the laboratory and finite‐difference data, we are able to propose a plausible explanation of the set of scattering events. Based on the laboratory data, we conclude that active seismic sources can potentially be used to determine fracture size and location in the field. The modeling example of fracture scattering illustrates the benefit of the finite‐difference technique with a variation in grid spacing for comparing numerical and physical experiments.

Optimal Decomposition of the Domain in Spectral Methods for Wave-Like Phenomena

A strategy for determining the optimal number of grid points and subdomains in a spectral method with domain decomposition on a serial computer is presented. The rapidly growing computational cost for large numbers of grid points in each subdomain is balanced against the exponential convergence for spectral approximation of smooth functions, and the optimum is found as the number of grid points and subdomains that gives the minimal computational cost for a given accuracy. The typical length scale of the problem is found to influence the number of subdomains but not the number of grid points within each subdomain.

Accurate viscoelastic modeling by frequency‐domain finite differences using rotated operators

The viscoelastic wave equation is an integro‐differential equation that requires special methods when using time‐domain numerical finite‐difference methods. In the frequency domain, the integral terms are easily represented by complex valued elastic media properties. There are further significant advantages to using the frequency domain if the forward or the inverse problem requires modeling or inverting a large number of prestack source gathers. Numerical modeling is expensive for seismic data because of the large number of wavelenghths typically separating sources from receivers, which results in a need for a large number of grid points. A major obstacle to using frequency‐domain methods is the consequent storage requirements. To reduce these, we maximize the accuracy and simultaneously minimize the spatial extent of the numerical operators. We achieve this by extending earlier published methods introduced for the viscoacoustic case to the viscoelastic case. This requires the formulation of two new numerical operators: a differencing operator in a rotated coordinate frame and a lumped mass term. The new operators are combined with ordinary second‐ order, finite‐difference operators in an optimal manner to minimize numerical errors without increasing the size of the numerical operator. For a fixed number of grid points, the resulting second‐ order differencing scheme is no more expensive than an ordinary second‐order differencing scheme, but a numerical dispersion analysis shows that the number of grid points required per smallest wavelength is reduced from approximately 15 to approximately 4. The new scheme is also capable of handling embedded fluid layers without instability. We demonstrate that no further improvement in performance can be achieved using higher order spatial operators because of the associated computational overheads associated with the larger differencing operators. The new viscoelastic modeling scheme is used to study a crosshole data set in which the exact nature of the seismic coda is unclear. The results of the modeling study indicate this coda is likely related to the generation of mode‐converted shear waves within the complicated, finely layered sediments at the site.

An Iterative PPPM Method for Simulating Coulombic Systems on Distributed Memory Parallel Computers

We describe results obtained from a new implementation of Hockney's Particle-Particle Particle-Mesh (PPPM) method for evaluation of Coulomb energies and forces in simulations of charged particles. Rather than taking the usual approach, solving Poisson's equation by means of a Fourier transformation, we use an iterative Poisson solver. In a molecular dynamics (MD) simulation the solution from the previous time-step provides a good starting point for the next solution. This reduces the number of iterations per time-step to acceptable values. The iterative scheme has a complexity O(N), and, in contrast with the Fourier transform based approach, it is easily implemented on a parallel architecture with a minimum of communication overhead. We examine the origin of the errors in the algorithm and find that reasonable accuracies in the Coulomb interaction can best be attained by making the charge density profile as smooth as possible. This involves spreading the particle charges over a large number of grid points. Assigning these charges then becomes the most time consuming part of the algorithm. We show how we can then gain a considerable saving in computing time by employing a diffusion equation as a charge spreading mechanism. The effect of employing the algorithm with an accuracy less than that typically tolerated in an Ewald summation is studied by computing, from an MD simulation of silica, quantities that are sensitive to the long range part of the Coulomb interaction. These results are compared to full Ewald sum reference simulations and found to be within the statistical error.