Hi,
I noticed a bug in the direct space part of ewald computation in
amber/sander (4.1 may be higher).
Direct calls to erfc are replaced by a cubic spline approximation:
c cubic spline on erfc,derf
ind = erftbdns*x + 1
dx = x - (ind-1)*del
erfcc = erf_arr(1,ind)+dx*(erf_arr(2,ind)+
$ dx*(erf_arr(3,ind) + dx*erf_arr(4,ind)*ONETHIRD)*0.5)
derfc = -(erf_arr(2,ind)+dx*(erf_arr(3,ind)+
$ dx*erf_arr(4,ind)* ONETHIRD))
With derfc, the first derivative of erfc is calculated by
differentiating the cubic spline by dx. To see it a bit more clear
f(x) = a1 + a2*(x) + 1/2*a3*(x^2) + 1/6*a4*(x^3)
f'(x)= a2*(1) + 1/2*a3*(2*x) + 1/6*a4*(3*x^2)
= a2 + a3*(x) + 1/2*a4*(x^2)
^^^
ONETHIRD in derfc is wrong
it should be 0.5 !!!
So every expression like
derfc = -(erf_arr(2,ind)+dx*(erf_arr(3,ind)+
$ dx*erf_arr(4,ind)* ONETHIRD))
has to be replaced by
derfc = -(erf_arr(2,ind)+dx*(erf_arr(3,ind)+
$ dx*erf_arr(4,ind)* 0.5d0))
Those changes might apply to ew_direct in sander and gibbs.
Is this bug corrected in amber 5.0?
Comments?
-----------------------------------------------------------------------------
Dr. Thomas Huber Institut fuer Stoffwechselbiochemie
Tel.: +49(89)5996-462 der Universitaet Muenchen
FAX: +49(89)5996-415 Vorstand Prof.Ch.Haass
email thuber_at_physik.tu-muenchen.de Schillerstrasse 44
D-80336 Muenchen
On Thu, 19 Aug 1999, Claudio Chuaqui wrote:
> Hi,
>
> In an effort to further investigate the Ewald energy differences
> on Intel processors from those obtained on SGI/IBM/SUN
> machines, I (with the help of others) found a bug in the gnu
> math library erfc function: /lib/libm-2.1.1.so.
>
> This has been submitted as a bug report to gnu and
> is now listed as closed on their database. The bug is not there
> in earlier versions of the library.
>
> The bug in the s_erf.c file occurs twice, once for erf and once for
> erfc. The top of the file reads,
>
> /* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/25,
> for performance improvement on pipelined processors.
> */
>
> so perhaps the error was introduced here. In any case, it seems that
> a coefficient, s6, is multiplied twice by error in 4 places. This
> is the diff (diff original fixed).
>
> ------------------------------------------
> diff /usr/src/redhat/SOURCES/glibc/sysdeps/libm-ieee754/s_erf.c.ORIG
> s_erf.c
>
> 249,250c249,250
> < P4 = s6*pa[6];
> < Q4 = s6*qa[6];
> ---
> > P4 = pa[6];
> > Q4 = qa[6];
> 365,366c365,366
> < P4 = s6*pa[6];
> < Q4 = s6*qa[6];
> ---
> > P4 = pa[6];
> > Q4 = qa[6];
> ------------------------------------------
>
> The bug can be verified easily by comparing erf and erfc values
> from the gnu library, at x=1.2 (for example), with those from
> the libraries on other machines. For example:
>
>
> derf derfc
> PIII 0.9103141590167326 8.9685840983267493E-002
>
> sun ultra 0.91031397822964 8.9686021770365D-02
> ibm 0.910313978229635334 8.96860217703646379E-02
> sgi 0.9103139782296354 8.9686021770364610E-02
>
> fixed 0.910313978229635334 8.9686021770364652
>
> With the error removed (fixed, above), the Pentium value is in
> agreement.
>
> Fixing this problem improved my Ewald energies such that they
> now agree much better with the test cases provided with
> Amber.
>
> Claudio Chuaqui
> Rowland Institute for Science
> Cambridge, MA
>
>
>
>
Received on Wed Aug 25 1999 - 05:52:21 PDT
