[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[ccp4bb]: Reading scattering factors (ATOMSF): bug in rwbrook.f



Dear all,

there is a bug in the CCP4 library (version 4.0.1 and probably earlier
ones too) for reading of atomic scattering factors (from $ATOMSF). It
occurs when calling subroutine SFREAD2 in $CCP4/lib/src/rwbrook.f to
get the 5 Guassian terms.

The problem occurs ONLY for 5 term Gaussian and ONLY for atoms S
(sulfur), I (iodine) and B (Boron). The atomic scattering factors
actually extracted for these atoms are:

  S  -->  Si
  I  -->  In
  B  -->  Be

(as far as I can tell, all other atoms should be fine). The 2 term
Gaussians for H, C, N, O and S are NOT affected. I think the following
programs are affected:

  AMORE
  MLPHARE
  SFALL
  TRUNCATE
  REFMAC (don't know about the latest version though ...)

I attach a quick fix for $CCP4/lib/src/rwbrook.f. The problem is, that
e.g. Si is before S in the periodic table (and also in $ATOMSF) and
therefore a single character match is found before reaching the proper
entry for S. The fix adds some special checks for single character
atom names.

Cheers

Clemens

PS: I also attach a patch to make SFALL print the sum[a(i)] in the
table of scattering factors - which should be the number of
electrons. That's an easy way to check if you're paranoid like me ...

-- 

***************************************************************
* Clemens Vonrhein, Ph.D.          vonrhein@GlobalPhasing.com
*
*  Global Phasing Ltd.
*  Sheraton House, Castle Park     Tel: +44-(0)1223-353033
*  Cambridge CB3 0AX, UK           Fax: +44-(0)1223-366889
*--------------------------------------------------------------
* BUSTER Development Group (http://Lagrange.mrc-lmb.cam.ac.uk)
***************************************************************
*** rwbrook.f.orig	Mon Nov 20 16:03:56 2000
--- rwbrook.f.new	Mon Nov 20 16:04:39 2000
***************
*** 2728,2733 ****
--- 2728,2737 ----
  C                 =  0 OK
  C                 =  1 for two gaussian case that does not exist
  C
+ C 20/11/2000 C. Vonrhein
+ C     fixed problem with one character atom name matching
+ C     (S, I and B affected)
+ C
  C     .. Scalar Arguments ..
        REAL C
        INTEGER IELEC,Ifail,IWT,NG
***************
*** 2781,2788 ****
  C
        CALL CCPUPC(IDIN)
        IF (ID2(1:NID).EQ.IDIN(1:NID)) THEN
!         Ifail = 1
!         IF (NGauss.NE.2 .OR. IDIN(6:6).NE.' ') GO TO 60
        END IF
  C
        GO TO 10
--- 2785,2804 ----
  C
        CALL CCPUPC(IDIN)
        IF (ID2(1:NID).EQ.IDIN(1:NID)) THEN
! c
! c       20/11/2000 C. Vonrhein
! c
! c       special precautions for single character atom types:
! c       skip this atom if second character is NOT '+', '-' or ' '.
! c
!         IF ((NID.NE.1).OR.( (NID.EQ.1).AND.(
!      .                      (IDIN(2:2).EQ."+").OR.
!      .                      (IDIN(2:2).EQ."-").OR.
!      .                      (IDIN(2:2).EQ." ")    )
!      .                    )) THEN
!           Ifail = 1
!           IF (NGauss.NE.2 .OR. IDIN(6:6).NE.' ') GO TO 60
!         END IF
        END IF
  C
        GO TO 10
*** /disk1/software/xpgm/src/ccp4_4.0/ccp4/src/dist/sfall.f	Tue Jan 18 10:39:27 2000
--- sfall_test.f	Mon Nov 20 15:45:18 2000
***************
*** 692,698 ****
        REAL AtmBoxLowLim,AtmBoxUpLim,CellMtz,ExtBoxLowLim,ExtBoxUpLim,
       +     ExtraPlanes,FormFactAcoeff,FormFactBcoeff,PermutRecipSymm,
       +     PlanesLimits,RealSymmMatrx,SpecialPlaneLimit,
!      +     RecipSymmMatrx
        INTEGER Iuvw,IuvwMP,KatmWghts,LSymmFlags,NumFirstResChn,
       +        NumGaussPerAtmTyp,NumLastResChn,NumResPerChn,Nxyz
        CHARACTER ChainLabels*1,ElementNames*4
--- 692,698 ----
        REAL AtmBoxLowLim,AtmBoxUpLim,CellMtz,ExtBoxLowLim,ExtBoxUpLim,
       +     ExtraPlanes,FormFactAcoeff,FormFactBcoeff,PermutRecipSymm,
       +     PlanesLimits,RealSymmMatrx,SpecialPlaneLimit,
!      +     RecipSymmMatrx,FormFactAcoeffSum
        INTEGER Iuvw,IuvwMP,KatmWghts,LSymmFlags,NumFirstResChn,
       +        NumGaussPerAtmTyp,NumLastResChn,NumResPerChn,Nxyz
        CHARACTER ChainLabels*1,ElementNames*4
***************
*** 747,753 ****
       +                NumGaussPerAtmTyp(MaxFormFactors),
       +                KatmWghts(MaxFormFactors),
       +                FormFactAcoeff(5,MaxFormFactors),
!      +                FormFactBcoeff(5,MaxFormFactors)
        COMMON /FORM2/
       +                ElementNames(MaxFormFactors)
        COMMON /HKLLIM/
--- 747,754 ----
       +                NumGaussPerAtmTyp(MaxFormFactors),
       +                KatmWghts(MaxFormFactors),
       +                FormFactAcoeff(5,MaxFormFactors),
!      +                FormFactBcoeff(5,MaxFormFactors),
!      +                FormFactAcoeffSum(MaxFormFactors)
        COMMON /FORM2/
       +                ElementNames(MaxFormFactors)
        COMMON /HKLLIM/
***************
*** 845,851 ****
            WRITE (6,FMT=6002)
   6002     FORMAT (///,
       +' Form factors used are   ',/,
!      +' AtomName - formfac number -  (AI BI),I=1,NGauss   ',/,
       +'  Default value is 2 term gaussian')
  C
            DO 10 N = 1,NumFormFactors
--- 846,852 ----
            WRITE (6,FMT=6002)
   6002     FORMAT (///,
       +' Form factors used are   ',/,
!      +' AtomName - formfac number -  (AI BI),I=1,NGauss - sum(AI)',/,
       +'  Default value is 2 term gaussian')
  C
            DO 10 N = 1,NumFormFactors
***************
*** 853,860 ****
       +                           KatmWghts(N),
       +                            (FormFactAcoeff(J,N),
       +                              FormFactBcoeff(J,N),
!      +                               J=1,NumGaussPerAtmTyp(N))
!  6004       FORMAT (2X,A4,I5,2X,2 (2F10.5,2X),/,9X,10 (F10.5,2X))
     10     CONTINUE
  C
            WRITE (6,FMT=6006) Title(1:LENSTR(Title))
--- 854,862 ----
       +                           KatmWghts(N),
       +                            (FormFactAcoeff(J,N),
       +                              FormFactBcoeff(J,N),
!      +                               J=1,NumGaussPerAtmTyp(N)),
!      +                                 FormFactAcoeffSum(N)
!  6004       FORMAT (2X,A4,I5,2X,2 (2F10.5,2X),/,13X,5 (2F10.5,2X))
     10     CONTINUE
  C
            WRITE (6,FMT=6006) Title(1:LENSTR(Title))
***************
*** 1376,1382 ****
        REAL AtmBoxLowLim,AtmBoxUpLim,CellMtz,ExtBoxLowLim,ExtBoxUpLim,
       +     ExtraPlanes,FormFactAcoeff,FormFactBcoeff,PermutRecipSymm,
       +     PlanesLimits,RealSymmMatrx,SpecialPlaneLimit,
!      +     RecipSymmMatrx
        INTEGER Iuvw,IuvwMP,KatmWghts,LookUp,LSymmFlags,NumFirstResChn,
       +        NumGaussPerAtmTyp,NumLastResChn,NumResPerChn,Nxyz
        CHARACTER ChainLabels*1,ElementNames*4
--- 1378,1384 ----
        REAL AtmBoxLowLim,AtmBoxUpLim,CellMtz,ExtBoxLowLim,ExtBoxUpLim,
       +     ExtraPlanes,FormFactAcoeff,FormFactBcoeff,PermutRecipSymm,
       +     PlanesLimits,RealSymmMatrx,SpecialPlaneLimit,
!      +     RecipSymmMatrx,FormFactAcoeffSum
        INTEGER Iuvw,IuvwMP,KatmWghts,LookUp,LSymmFlags,NumFirstResChn,
       +        NumGaussPerAtmTyp,NumLastResChn,NumResPerChn,Nxyz
        CHARACTER ChainLabels*1,ElementNames*4
***************
*** 1463,1469 ****
       +                NumGaussPerAtmTyp(MaxFormFactors),
       +                KatmWghts(MaxFormFactors),
       +                FormFactAcoeff(5,MaxFormFactors),
!      +                FormFactBcoeff(5,MaxFormFactors)
        COMMON /FORM2/
       +                ElementNames(MaxFormFactors)
        COMMON /IOSF/
--- 1465,1472 ----
       +                NumGaussPerAtmTyp(MaxFormFactors),
       +                KatmWghts(MaxFormFactors),
       +                FormFactAcoeff(5,MaxFormFactors),
!      +                FormFactBcoeff(5,MaxFormFactors),
!      +                FormFactAcoeffSum(MaxFormFactors)
        COMMON /FORM2/
       +                ElementNames(MaxFormFactors)
        COMMON /IOSF/
***************
*** 1640,1645 ****
--- 1643,1650 ----
          FormFactBcoeff(3,IFMFC) = BF(3)
          FormFactBcoeff(4,IFMFC) = BF(4)
          FormFactAcoeff(5,IFMFC) = CF
+         FormFactAcoeffSum(IFMFC)=
+      +    AF(1) + AF(2) + AF(3) + AF(4) + CF
          ElementNames(IFMFC)   = ATID(1:4)
    5   CONTINUE
  C
***************
*** 2374,2379 ****
--- 2379,2386 ----
          FormFactBcoeff(3,NumFormFactors) = BF(3)
          FormFactBcoeff(4,NumFormFactors) = BF(4)
          FormFactAcoeff(5,NumFormFactors) = CF
+         FormFactAcoeffSum(NumFormFactors)= 
+      +    AF(1) + AF(2) + AF(3) + AF(4) + CF
          ElementNames(NumFormFactors)   = ATID(1:4)
    315 CONTINUE
        END IF
***************
*** 2419,2424 ****
--- 2426,2433 ----
          FormFactBcoeff(3,NumFormFactors)    = BF(3)
          FormFactBcoeff(4,NumFormFactors)    = BF(4)
          FormFactAcoeff(5,NumFormFactors)    = CF
+         FormFactAcoeffSum(NumFormFactors)   = 
+      +    AF(1) + AF(2) + AF(3) + AF(4) + CF
          ElementNames(NumFormFactors)      = ATID(1:4)
    320 CONTINUE
  C
***************
*** 3771,3777 ****
        REAL AtmBoxLowLim,AtmBoxUpLim,CellMtz,ExtBoxLowLim,ExtBoxUpLim,
       +     ExtraPlanes,PermutRecipSymm,PlanesLimits,
       +     FormFactAcoeff,FormFactBcoeff,
!      +     RealSymmMatrx,SpecialPlaneLimit,RecipSymmMatrx
        INTEGER Iuvw,LSymmFlags,NumFirstResChn,NumLastResChn,NumResPerChn,
       +     KatmWghts,NumGaussPerAtmTyp,Ifail,iwt,ielec
        CHARACTER ChainLabels*1,ElementNames*4,MAACD1*1,MAACID*4,MATMTP*4
--- 3780,3787 ----
        REAL AtmBoxLowLim,AtmBoxUpLim,CellMtz,ExtBoxLowLim,ExtBoxUpLim,
       +     ExtraPlanes,PermutRecipSymm,PlanesLimits,
       +     FormFactAcoeff,FormFactBcoeff,
!      +     RealSymmMatrx,SpecialPlaneLimit,RecipSymmMatrx,
!      +     FormFactAcoeffSum
        INTEGER Iuvw,LSymmFlags,NumFirstResChn,NumLastResChn,NumResPerChn,
       +     KatmWghts,NumGaussPerAtmTyp,Ifail,iwt,ielec
        CHARACTER ChainLabels*1,ElementNames*4,MAACD1*1,MAACID*4,MATMTP*4
***************
*** 3831,3837 ****
       +                NumGaussPerAtmTyp(MaxFormFactors),
       +                KatmWghts(MaxFormFactors),
       +                FormFactAcoeff(5,MaxFormFactors),
!      +                FormFactBcoeff(5,MaxFormFactors)
        COMMON /FORM2/
       +                ElementNames(MaxFormFactors)
        COMMON /IOSF/
--- 3841,3848 ----
       +                NumGaussPerAtmTyp(MaxFormFactors),
       +                KatmWghts(MaxFormFactors),
       +                FormFactAcoeff(5,MaxFormFactors),
!      +                FormFactBcoeff(5,MaxFormFactors),
!      +                FormFactAcoeffSum(MaxFormFactors)
        COMMON /FORM2/
       +                ElementNames(MaxFormFactors)
        COMMON /IOSF/
***************
*** 4092,4108 ****
          FormFactBcoeff(3,NumFormFactors)    = BF(3)
          FormFactBcoeff(4,NumFormFactors)    = BF(4)
          FormFactAcoeff(5,NumFormFactors)    = CF
          ElementNames(NumFormFactors)      = ATID
            WRITE (6,'(/,I4,a,/,a)')
       +NumFormFactors,'th Form factor added:   ',
!      +' AtomName - atomic number  -  (AI BI),I=1,5   '
  C
!           WRITE (6,'(2X,A4,I15,2X,2 (2F10.5,2X),/,11X,3 (2F10.5,2X),/)')
!      +                          ElementNames(NumFormFactors),
!      +                           KatmWghts(NumFormFactors),
!      +                            (FormFactAcoeff(J,NumFormFactors),
!      +                              FormFactBcoeff(J,NumFormFactors),
!      +                               J=1,5)
     26       CONTINUE
  c          END IF 
  
--- 4103,4122 ----
          FormFactBcoeff(3,NumFormFactors)    = BF(3)
          FormFactBcoeff(4,NumFormFactors)    = BF(4)
          FormFactAcoeff(5,NumFormFactors)    = CF
+         FormFactAcoeffSum(NumFormFactors)   = 
+      +    AF(1) + AF(2) + AF(3) + AF(4) + CF
          ElementNames(NumFormFactors)      = ATID
            WRITE (6,'(/,I4,a,/,a)')
       +NumFormFactors,'th Form factor added:   ',
!      +' AtomName - atomic number  -  (AI BI),I=1,5  -  sum(AI) '
  C
!           WRITE (6,6005) ElementNames(NumFormFactors),
!      +                    KatmWghts(NumFormFactors),
!      +                     (FormFactAcoeff(J,NumFormFactors),
!      +                       FormFactBcoeff(J,NumFormFactors),
!      +                        J=1,5),
!      +                         FormFactAcoeffSum(NumFormFactors)
!  6005     FORMAT(2X,A4,I15,2X,2 (2F10.5,2X),/,23X,3(2F10.5,2X),F10.5,/)
     26       CONTINUE
  c          END IF 
  
***************
*** 7094,7100 ****
  C     .. Arrays in Common ..
        REAL AtmBoxLowLim,AtmBoxUpLim,CellMtz,ExtBoxLowLim,ExtBoxUpLim,
       +     ExtraPlanes,FormFactAcoeff,FormFactBcoeff,PermutRecipSymm,
!      +     PlanesLimits,RealSymmMatrx,SpecialPlaneLimit,RecipSymmMatrx
        INTEGER Iuvw,KatmWghts,LSymmFlags,NumGaussPerAtmTyp
        CHARACTER MAACD1*1,MAACID*4,MATMTP*4
  C     ..
--- 7108,7115 ----
  C     .. Arrays in Common ..
        REAL AtmBoxLowLim,AtmBoxUpLim,CellMtz,ExtBoxLowLim,ExtBoxUpLim,
       +     ExtraPlanes,FormFactAcoeff,FormFactBcoeff,PermutRecipSymm,
!      +     PlanesLimits,RealSymmMatrx,SpecialPlaneLimit,RecipSymmMatrx,
!      +     FormFactAcoeffSum
        INTEGER Iuvw,KatmWghts,LSymmFlags,NumGaussPerAtmTyp
        CHARACTER MAACD1*1,MAACID*4,MATMTP*4
  C     ..
***************
*** 7149,7155 ****
       +                NumGaussPerAtmTyp(MaxFormFactors),
       +                KatmWghts(MaxFormFactors),
       +                FormFactAcoeff(5,MaxFormFactors),
!      +                FormFactBcoeff(5,MaxFormFactors)
        COMMON /IOSF/
       +                IOscratch,IOgradmat,IOhessian,IOshifts,
       +                IOatomop,IOmap,IOxyzin,IOxyzoutunq
--- 7164,7171 ----
       +                NumGaussPerAtmTyp(MaxFormFactors),
       +                KatmWghts(MaxFormFactors),
       +                FormFactAcoeff(5,MaxFormFactors),
!      +                FormFactBcoeff(5,MaxFormFactors),
!      +                FormFactAcoeffSum(MaxFormFactors)
        COMMON /IOSF/
       +                IOscratch,IOgradmat,IOhessian,IOshifts,
       +                IOatomop,IOmap,IOxyzin,IOxyzoutunq
***************
*** 7692,7698 ****
        LOGICAL LNoScaleFlag,LverboseFlag
  C     ..
  C     .. Arrays in Common ..
!       REAL CellMtz,FormFactAcoeff,FormFactBcoeff
        INTEGER Iuvw,KatmWghts,NumGaussPerAtmTyp
  C     ..
  C     .. Local Scalars ..
--- 7708,7714 ----
        LOGICAL LNoScaleFlag,LverboseFlag
  C     ..
  C     .. Arrays in Common ..
!       REAL CellMtz,FormFactAcoeff,FormFactBcoeff,FormFactAcoeffSum
        INTEGER Iuvw,KatmWghts,NumGaussPerAtmTyp
  C     ..
  C     .. Local Scalars ..
***************
*** 7720,7726 ****
       +                NumGaussPerAtmTyp(MaxFormFactors),
       +                KatmWghts(MaxFormFactors),
       +                FormFactAcoeff(5,MaxFormFactors),
!      +                FormFactBcoeff(5,MaxFormFactors)
        COMMON /FORM2/
       +                ElementNames(MaxFormFactors)
        COMMON /HKLLIM/
--- 7736,7743 ----
       +                NumGaussPerAtmTyp(MaxFormFactors),
       +                KatmWghts(MaxFormFactors),
       +                FormFactAcoeff(5,MaxFormFactors),
!      +                FormFactBcoeff(5,MaxFormFactors),
!      +                FormFactAcoeffSum(MaxFormFactors)
        COMMON /FORM2/
       +                ElementNames(MaxFormFactors)
        COMMON /HKLLIM/