[ccp4bb]: Bug in numerical support library for refmac

[UW/BMSC CCP4 users <ccp4@gouda.bmsc.washington.edu>]

***  For details on how to be removed from this list visit the  ***
***    CCP4 home page http://www.dl.ac.uk/CCP/CCP4/main.html    ***

Dear fellow CCP4 refiners:

I have found a bug that causes refmac, both versions 4 and 5, to crash
on a divide-by-zero error under some circumstances.  This bug is only
triggered during refinement of individual anisotropic parameters Uij,
and the specific input values that trigger it probably vary with the
computer architecture and compiler options.  Nevertheless it is a bug,
and I was only able to proceed with my current aniso refinement after
tracking it down and squashing it.  

I recommend applying the fix below to all versions of refmac, and also 
to anisoanal.f, molrep_prog.f, and anywhere else that the EIGEN_ASC
routine appears in the current CCP4 distribution.

In routine EIGEN_ASC of aniso_refi.f the following lines appear 
(line 2084 of aniso_refi.f in version 4;
 line 1851 of aniso_refi_hd1.f in version 5)
	C
	C---- COMPUTE SIN & COS
	60    MQ=M*(M-1)/2
	      LM=L+MQ
        CEAM  IF (ABS(A(LM))-THR.LT.0.) GOTO 130
	      IF (A(LM)*A(LM)-THR.LT.0.) GOTO 130
	      IND=1
	      MM=M+MQ
	      X=.5*(A(LL)-A(MM))
	      Y=-A(LM)/SQRT(A(LM)**2+X*X)

Note that I have commented out a test on abs(A(LM)) and replaced it
with a test on A(LM)**2.  This is necessary because of the specific
case where X=0 and A(LM) is sufficiently small that A(LM)**2 causes
floating underflow and gets rounded to 0.0.  In that case the line
Y = -A(LM) / sqrt(A(LM)**2 + X*X)  triggers a divide-by-zero error
and the program crashes.  

The deeper cause of the problem is that the threshold value for
numerical convergence, THR, is becoming much too small.
But fixing that would require more extensive code changes and I didn't
want to fiddle with what looks to be ancient IBM library code.
If improved performance and general architecture independence is an 
issue here, it would probably be better to replace all these support 
routines with calls to the equivalent routines in the lapack library, 
which is actively maintained and is tweaked during installation to 
match the host architecture and compiler.

			regards,

				Ethan A Merritt

-------------------------------------------------------------------
Ethan A Merritt                         K428b Health Sciences
Dept of Biological Structure            phone: (206)543-1421
Mailstop 357742         		FAX:   (206)685-7002
University of Washington
Seattle, WA 98195-7742                  merritt@u.washington.edu
-------------------------------------------------------------------