# [gmx-users] PME accuracy

Vincent Ballenegger vcb25 at cam.ac.uk
Thu May 8 17:55:01 CEST 2003

```Erik Lindahl wrote:

>First - the accuracy estimates are those of Darden and coworkers; check
>out the Smooth PME paper mentioned in the references.
>
>Second, I'll have a look both at PME and the normal (non-PME) fourier
>method. We haven't used that for 1-2 years, so there might be a bug
>there rather than in PME...
>
>
Dear all,
I had a closer look at the smooth PME method, and I think my results
might interest all users of the PME method.

I compared the results of the classical Ewald method and the smooth PME
method for the electrostatic force between two ions (charge +e and -e)
along the diagonal of a cubix box of length 3.5nm, for interionic
distances between 0.35 and 3mn.
I measured three curves (see plot), with the following parameters:
- Ewald method: rcoulomb = 0.9, ewald_rtol = 1e-5, fourier_nx,y,z = 25
- PME method: rcoulomb = 0.9, ewald_rtol = 1e-5, fourier_spacing = 0.1
(=> grid 35 x 35 x 35)
pme_order = 4 for the first curve
pme_order = 6 for the second curve
We see immediately on the plot that
1) the classical Ewald method works and gives a nice smooth curve.
(I checked that this curve converges to the Coulomb force 1/r^2 when the
box size increases.)
2) the PME curve with pme_order = 6 falls nicely on the previous curve
3) the PME curve with pme_order = 4 is reasonably accurate at short
distances, but develops rapidely large oscillations and gets quickly
very inaccurate. For example, the relative error in the force is 6% at
1.5nm, 42% at 2.5nm, 110% at 2.84nm !!

I checked the paper on the smooth PME method, and their estimation of
its accuracy (they claim a RMS error of 4e-4 with pme_order=4 in their
Table I). I think their estimate of the accuracy is very misleading
because they don't calculate the true RMS error in the full (i.e.
direct+reciprocal) electrostatic forces. Instead of defining

RMS_error = sqrt{ Sum[(F_pme - F_ewald)^2] } /  average(F_ewald)

they use

reciprocal_error = sqrt{ Sum[(F_pme_reciprocal -
F_ewald_reciprocal)^2]  /  sum(F_ewald^2) }.

Yes, they really do compute the sum of the squares in the denominator,
instead of the average of F_ewald, and then take the square root!!! (For
some strange reason, they use also F_ewald^2 instead of
F_ewald_reciprocal^2 in the denominator.)

Since F_ewald is very large at short distances, it is clear that
reciprocal_error is a small number (about 1e-3). The true RMS error is
however much larger ! The RMS errors for two PME curves I calculated are
as follows:

pme_order 4                 pme_order 6

a) RMS_error = 0.7%         a) RMS_error = 0.02%
b) RMS_error = 6.4%         b) RMS_error = 0.21%
c) RMS_error = 43.1%        c) RMS_error = 1.45%

The average (a) is calculated for distances between 0 and 1.17nm (first
third of the curve).
The average (b) is for the second third of the curve (1.17<d<2.33nm),
and (c) is for d>2.33nm.

In conclusion:  SETTING pme_order TO 4 LEADS TO VERY LARGE ERRORS IN THE
FORCES.

Gromacs default value for pme_order should therefore be at least 6