[gmx-users] PME accuracy

[Vincent Ballenegger]

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
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: PLOT.agr
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20030508/aa1ef5a1/attachment.ksh>