Dear Thomas,
You are right about the bug. I think it has a small effect; usually in the
6th significant figure, which is why it escaped notice until recently. It
was fixed some time between amber5 and 6. Its there in my amber5.0 source
but not in amber6 (soon to come out)
Tom Darden
On Wed, 25 Aug 1999, Thomas Huber wrote:
> 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 - 06:58:09 PDT
