[gmx-users] PME accuracy

[Jason DeJoannis]

Vincent,

  The error estimate of 4e-4 is essentially that: an estimate 
-- and a crude one at that. A thorough analysis of errors 
was done by Kolafa & Perram (Mol Sim, 1992) and by Petersen 
(JCP, 1995). They find that the error does not depend too 
much on the configuration for __large__ N. As you have 
noticed it does depend on the configuration for small N. 
Don't worry, the long-distance, large relative-errors 
you noted will 'wash out' when the number of charges is 
increased, as there are necessarily other particles nearby 
that make a bigger contribution to the force on the particle 
of interest. The presence of a cosine is no surprise either. 
Look at the regular Ewald equation for the Fourier space. 
There is a cosine there, and it will show up under special 
circumstances. See for example the two-particle case I showed 
in JCP, v. 117, p. 2503, fig. 3. We too were puzzled by this
oddity for a while.

  As a rule we should keep in mind that 4e-4 is a relative-error 
estimate for "typical" systems. That is, systems with a 
typical size (~4nm) and a typical number of charges (~6000)
with a non-problematic spatial distribution. In fact I cannot
think of any problematic distributions so the latter condition
is really not of concern.

  For more of the nitty-gritty read on:

  The distribution of errors in a many-particle system has
been studied. Petersen even defined a worst-case error and
this is well-behaved. It should be well-behaved because the
truncation of the Ewald sum for a given particle produces
a Gaussian error distribution (again see the above 3 refs). 
I haven't heard any reason why the spline or polynomial 
interpolation should change this fact.

  There are a few situations when the the worst-case becomes
an issue (i.e. much larger than RMS). To my knowledge this
only happens for 2D methods: MMM2D, EW3DC and ELC and for 
particles near the wall. For these situations "maximal error" 
formulas were derived (fairly rigorously).

  The best way to check that your installation of PME is
working properly is probably to use a box of water molecules.
As you have done, one can compare the results to a well-converged
Ewald sum. You should get roughly 4e-4 for any water 
configuration. If you are in the mood for a more exact 
test (and more work) you can use the 100-particle system 
of Deserno & Holm (JCP, 1998) and compare to many
decimal-places.

Best regards,

  -Jason

---
Jason de Joannis, Ph.D.
Chemistry Department, Emory University
1515 Pierce Dr. NE, Atlanta, GA 30322
Phone: (404) 712-2983
Email: jdejoan at emory.edu
http://userwww.service.emory.edu/~jdejoan



Quoting Vincent Ballenegger <vcb25 at cam.ac.uk>:

> 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 
> (leading to RMS errors of about 1%).
> 
> Cheers,
> 
> Dr. Vincent Ballenegger
> University of Cambridge
>