[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
(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>
More information about the gromacs.org_gmx-users
mailing list