Re: [ccp4bb]: SFALL and anomalous dispersion

```Somewhat laborious - we can calculat Fano from the anomalous scatterers
then add together the real and inmaginary SFS with appropriate
multiplication factors

eg F(+) = Ftot -f'*Fano +iF" Fano
F(-) = Ftot* -f'Fano* +if"Fano*  where F* is the complex conjugate

It uses a jiffy called mtzadd.f
I will attach scripts - all pre GUI days!

Elranor
And mtzadd.f - not sure if it is in distribution..

Bernard Santarsiero wrote:
> ***  For details on how to be removed from this list visit the  ***
>
> Is there a way to calculate F(+) and F(-) or DANO in SFALL or CPP4?
>
> Bernie Santarsiero
>
>

```
```#!/bin/csh -f
#
#   Calculate structure factors for all Atase atoms include Ses;
#   standard formfactors
goto 2
sfall		\
HKLIN /x/people/yao/atase/ecse3_sort_trun_l123sc.mtz \
HKLOUT /y/ccp4/work/1ecf_x0.5_4.mtz \
XYZIN ./1ecf_x0.5_4.pdb \
<< 'END-sfrkall'
TITL Structure factors calculed for all SD atoms constant formactor
MODE SFCALC XYZIN   HKLIN
RESO 37 2.25
BINS  60
LABI FP=FSe1 SIGFP=SIGFSe1
LABO ALLIN FC=FC_x0.5_4 PHIC=AC_x0.5_4
END
'END-sfrkall'
#
#
#  Get structure factors for the Se alone, using a Unitary structure factor
# Extract from Se pdb file
#  CRYST1  116.900  157.500  106.300  90.00  90.00  90.00 C 2 2 21     16
#  SCALE1      0.008554  0.000000  0.000000        0.00000
#  SCALE2      0.000000  0.006349  0.000000        0.00000
#  SCALE3      0.000000  0.000000  0.009407        0.00000
#  ATOM     67 An   MET A  10  1   27.688  55.343  19.523  1.00 27.19
#  ATOM    455 An   MET A  60  1   14.303  49.818  24.775  1.00 17.68

#  Extract from atomsf.lib for the unitary formfactor
#
#  An   2
#         102       102        0.000000
#          1.000000        0.000000        0.000000        0.000000
#          0.000000        0.000000        0.000000        0.000000
#          0.000000        1.000000        0.000000        1.000000

sfall		\
HKLIN  /y/ccp4/work/1ecf_x0.5_4.mtz \
HKLOUT /y/ccp4/work/1ecf_sd_x0.5.mtz \
XYZIN ./1ecf_sd_x0.5_4.pdb \
ATOMSF /y/ccp4/atase/atomsf.lib \
<< 'END-sfrkall'
TITL Structure factors calculed for all SD atoms constant formactor
MODE SFCALC XYZIN HKLIN
RESO 37 2.25
NOSCALE
BINS  60
LABI FP=FSe1 SIGFP=SIGFSe1
LABO ALLIN FC=FC_seu PHIC=AC_seu
END
'END-sfrkall'
#
#
#!/bin/csh -f
#
2:
hklin  /y/ccp4/work/1ecf_sd_x0.5.mtz \
hklout  /y/ccp4/work/junk1.mtz \
<<eof
LABI FC1=FC_x0.5_4 PHIC1=AC_x0.5_4 -
FC2=FC_seu    PHIC2=AC_seu -
FC3=FC_seu    PHIC3=AC_seu
LABO FSUM = FCL1_x0.5_4(+) PHISUM=PHICL1_x0.5_4(+)
SCALE  FC1 1 0                #   Scale Bfactor
SCALE  FC2 -5.5 0             #   Scale Bfactor
SCALE  COMP FC3 1.3 0         #   Complex ( ie advanvce phase by 90) Scale Bfactor
END
eof

#
hklin   /y/ccp4/work/junk1.mtz \
hklout  /y/ccp4/work/junk2.mtz \
<<eof
#LABI FC1=FC_x0.5_4 PHIC1=AC_x0.5_4 FC2=FC_seu PHIC2=AC_seu
LABI FC1=FC_x0.5_4 PHIC1=AC_x0.5_4 -
FC2=FC_seu    PHIC2=AC_seu -
FC3=FC_seu    PHIC3=AC_seu
LABO FSUM = FCL1_x0.5_4(-) PHISUM=PHICL1_x0.5_4(-)
SCALE  FC1 1 0
SCALE  FC2 -5.5 0
SCALE  COMP FC3 -1.3 0
END
eof

#
hklin   /y/ccp4/work/junk2.mtz \
hklout  /y/ccp4/work/junk3.mtz \
<<eof
LABI FC1=FC_x0.5_4 PHIC1=AC_x0.5_4 -
FC2=FC_seu    PHIC2=AC_seu -
FC3=FC_seu    PHIC3=AC_seu
LABO FSUM = FCL2_x0.5_4(+) PHISUM=PHICL2_x0.5_4(+)
SCALE  FC1 1 0
SCALE  FC2 -4.5 0
SCALE  COMP FC3 2.4 0
END
eof

#
hklin   /y/ccp4/work/junk3.mtz \
hklout  /y/ccp4/work/junk4.mtz \
<<eof
LABI FC1=FC_x0.5_4 PHIC1=AC_x0.5_4 -
FC2=FC_seu    PHIC2=AC_seu -
FC3=FC_seu    PHIC3=AC_seu
LABO FSUM = FCL2_x0.5_4(-) PHISUM=PHICL2_x0.5_4(-)
SCALE  FC1 1 0
SCALE  FC2 -4.5 0
SCALE  COMP FC3 -2.4 0
END
eof

#
hklin   /y/ccp4/work/junk4.mtz \
hklout  /y/ccp4/work/junk5.mtz \
<<eof
LABI FC1=FC_x0.5_4 PHIC1=AC_x0.5_4 -
FC2=FC_seu    PHIC2=AC_seu -
FC3=FC_seu    PHIC3=AC_seu
LABO FSUM = FCL3_x0.5_4(+) PHISUM=PHICL3_x0.5_4(+)
SCALE  FC1 1 0
SCALE  FC2 -2.0 0
SCALE  COMP FC3 1.4 0
END
eof

#
hklin   /y/ccp4/work/junk5.mtz \
hklout  /y/ccp4/work/junk6.mtz \
<<eof
LABI FC1=FC_x0.5_4 PHIC1=AC_x0.5_4 -
FC2=FC_seu    PHIC2=AC_seu -
FC3=FC_seu    PHIC3=AC_seu
LABO FSUM = FCL3_x0.5_4(-) PHISUM=PHICL3_x0.5_4(-)
SCALE  FC1 1 0
SCALE  FC2 -2.0 0
SCALE  COMP FC3 -1.4 0
END
eof
#
hklout ./1ecf_x0.5_4_L1L2L3.mtz \
<<eof
LABI -
F1+=FCL1_x0.5_4(+)  SIGF1+=SIGFSe(+)1 F1-=FCL1_x0.5_4(-) SIGF1-=SIGFSe(-)1 -
F2+=FCL2_x0.5_4(+)  SIGF2+=SIGFSe(+)2 F2-=FCL2_x0.5_4(-) SIGF2-=SIGFSe(-)2 -
F3+=FCL3_x0.5_4(+)  SIGF3+=SIGFSe(+)3 F3-=FCL3_x0.5_4(-) SIGF3-=SIGFSe(-)3

LABO -
F1=FCL1_x0.5_4 SIGF1=SIGFCL1_x0.5_4 D1=DCL1_x0.5_4 SIGD1=SIGDCL1_x0.5_4 -
F2=FCL2_x0.5_4 SIGF2=SIGFCL2_x0.5_4 D2=DCL2_x0.5_4 SIGD2=SIGDCL2_x0.5_4 -
F3=FCL3_x0.5_4 SIGF3=SIGFCL3_x0.5_4 D3=DCL3_x0.5_4 SIGD3=SIGDCL3_x0.5_4
END
eof
exit

```
```C
C     This code is distributed under the terms and conditions of the
C     CCP4 licence agreement as `Part ii)' software.  See the conditions
C     in the CCP4 manual for a copyright statement.
C
C     Authors:  EJD, MDW
C
C  mtzmnf.f
C STANDARD ITEM NAMES taken from \$CDOC/mtzlib.doc
C        Standard names are normally used for the following items:
C
C       Name          Item
C
C       H, K, L        Miller indices.
C   .....
C       I              Intensity.
C       SIGI           sigma(I).
C   .....
C       FP             Native 'F' value.
C       FC             Calculated 'F'.
C       FPHn           'F' value for derivative 'n'.
C
C       DP             Anomalous difference for native data.
C       DPHn           Anomalous difference for derivative 'n'.
C
C       SIGFP          sigma(FP).
C       SIGDP          sigma(DP).
C       SIGFPHn        sigma(Fn).
C       SIGDPHn        sigma(DELn).
C
C       PHIC           Calculated Phase.
C       PHIM           Most probable phase.
C       PHIB           Phase.
C
C       FOM            figure of merit.
C       WT             weight
C

CCOLUMN TYPES
C       The column types are as follows:
C
C       H       index h,k,l
C       J       intensity
C       F       structure amplitude, F
C       D       anomalous difference
C       Q       standard deviation of anything: J,F,D or other
C       P       phase angle in degrees
C       W       weight (of some sort )
C       A       phase probability coefficients (Hendrickson/Lattmann)
C   .....
C       I       any other integer
C       R       any other real

C*************************************************************************************

C---------  description of Key_Worded Input --------------
C           ================================
C
C
C  (i)   TITLE      (optional)
C
C        Replaces title in MTZ file
C
C  (ii)  LABIN      (compulsory)
C
C        Examples:
C           LABI F1=FTOXD3 SIGF1=SIGFTOXD3 -
C                F2=FAU20  SIGF2=SIGFAU20 -
C                D2=ANAU20 SIGD2=SIGANAU20 -
C
C  (iii) END      (optional)
C
C        Ends keyword input
C
PROGRAM MAIN
C     ..
C     .. Parameters ..
INTEGER MCOLS,NKEYWD,LSYM
PARAMETER (MCOLS=200,NKEYWD=200,LSYM=192)
C     ..
C     .. Local Arrays ..
+     SCL(MCOLS),B(MCOLS),COMPL(MCOLS)
INTEGER LOOKUP(MCOLS),IN(3),NCHANGES(MCOLS)
CHARACTER*1 CTPRGI(MCOLS),CTUSRI(MCOLS)*1,CTPRGO(2)
CHARACTER*30 LSPRGI(MCOLS),CLABS(MCOLS),LSUSRI(MCOLS)*30,
+             LSPRGO(2)
LOGICAL LOGMSS(MCOLS),LMSSOLD(MCOLS)
C     ..
C     .. Local Scalars ..
REAL S
INTEGER I,IC,IFAIL,IFIL1,IFLAG,IPRINT,JDO30,JDO35,JDO40,JDO60
INTEGER JDO65,JDO70,JDO80,JDO90,LUNOUT,NCOL,NF1,NLPRGI,NSPGRP
INTEGER NSYM,NSYMP,NUMFO,NUMFC,NUMI,NUMPH,NTOK,STARTFO,STARTFC
INTEGER STARTI,STARTPH
CHARACTER NAMSPG*10,PGNAM*10,TITLE*80
CHARACTER LTYPEX*1
LOGICAL ENDFIL,LTITLE
C     ..
C     .. Variables for Parser ..
REAL FVALUE(NKEYWD)
INTEGER IBEG(NKEYWD),IDEC(NKEYWD),IEND(NKEYWD),ITYPE(NKEYWD)
CHARACTER KEY*4,CVALUE(NKEYWD)*4,LINE*600,WORD*4
LOGICAL LEND
C     ..
C     .. External Routines ..
EXTERNAL CCPERR,CCPFYP,CCPRCS,CENTR,CENTRIC,EQUAL_MAGIC,LKYIN,
+         LROPEN,LRASSN,LRCLAB,LRCLOS,LRREFL,LRREFM,LRSYMI,LRSYMM,
+         LWOPEN,LWREFL,MTZINI,PARSER
C     ..
C     .. Data Statements ..
DATA NLPRGI/43/
DATA NLPRGO/2/
DATA LTITLE/.FALSE./
C  Allow NUMFO sets of Fs and phases
C  Up to NUMFC FCs
C  Up to NUMI Is
DATA NUMFO,NUMFC,NUMI/5,5,5/
DATA LSPRGI/'H','K','L',
+      'F1' ,'SIGF1' ,'PHIB1' ,'SIGD1' ,
+      'F2' ,'SIGF2' ,'PHIB2' ,'SIGD2' ,
+      'F3' ,'SIGF3' ,'PHIB3' ,'SIGD3' ,
+      'F4' ,'SIGF4' ,'PHIB4' ,'SIGD4' ,
+      'F5' ,'SIGF5' ,'PHIB5' ,'SIGD5' ,
+      'F6' ,'SIGF6' ,'PHIB6' ,'SIGD6' ,
+      'FC1','PHIC1' ,
+      'FC2','PHIC2' ,
+      'FC3','PHIC3' ,
+      'FC4','PHIC4' ,
+      'FC5','PHIC5' ,
+      'I1' ,'SIGI1' ,
+      'I2' ,'SIGI2' ,
+      'I3' ,'SIGI3' ,
+      'I4' ,'SIGI4' ,
+      'I5' ,'SIGI5' ,
+     153*' '/
DATA CTPRGI/'H','H','H',
+      'F','Q','P','W',
+      'F','Q','P','W',
+      'F','Q','P','W',
+      'F','Q','P','W',
+      'F','Q','P','W',
+      'F','Q','P','W',
+      'F','P',
+      'F','P',
+      'F','P',
+      'F','P',
+      'F','P',
+      'J','Q',
+      'J','Q',
+      'J','Q',
+      'J','Q',
+      'J','Q',
+     153*' '/
DATA LSPRGO/'FSUM','PHISUM'/
DATA CTPRGO/'F','P'/
DATA LOOKUP/-1,-1,-1,0,0,0,0,0,0,0,0,189*0/
DATA NCHANGES/200*0/
C
C     Starting indices of program label blocks
STARTFO=4
STARTFC=NUMFO*4 + STARTFO
STARTI=NUMFC*2 + STARTFC
LUNOUT=6
IFIL1 = 1
NREF = 0
PI=3.14159265
DTOR=PI/180.

CALL CCPFYP
CALL MTZINI

20   CONTINUE
KEY  = ' '
LINE = ' '
C     Reset number of words on line from 20 to NKEYWD
C     Need lots for LABI lines
NTOK = NKEYWD
C
C          **********************************************************
CALL PARSER(KEY,LINE,IBEG,IEND,ITYPE,FVALUE,CVALUE,IDEC,NTOK,LEND,
+     .TRUE.)
C          **********************************************************
C
IF (LEND) GO TO 200

IF(KEY.EQ.'TITL') THEN
IF (NTOK.GE.2) THEN
LTITLE = .TRUE.
TITLE = LINE(IBEG(2) :IEND(NTOK))
WRITE(LUNOUT,'(10X,a,/,15X,a80,/)')
+                        'TITLE OF MTZ FILE REPLACED BY:',TITLE
ENDIF
GO TO 20

ELSE IF(KEY.EQ.'LABI') THEN
C----    LABIN     i.e. set input labels
C        *****************************************
CALL LKYIN(IFIL1,LSPRGI,NLPRGI,NTOK,LINE,IBEG,IEND)
C        *****************************************
GO TO 20

ELSE IF(KEY.EQ.'LABO') THEN
C----    LABOUT    i.e. set output labels
C        *****************************************
CALL LKYOUT(IFIL1,LSPRGO,NLPRGO,NTOK,LINE,IBEG,IEND)
C        *****************************************
GO TO 20

C---- SCALE  [Complex] Label Scale B-factor
C     ===========================
C
C           F1 or F2 Scale and temperature factor
C
C
ELSE IF(KEY.EQ.'SCAL') THEN
ITOK = 2
COMPLEX = 0.0
WORD = LINE(IBEG(2):IBEG(2))
CALL CCPUPC(WORD)
IF(WORD(1:1) .EQ.'C' ) THEN
COMPLEX = 1.0
ITOK = 3
END IF
5        IF (ITOK .LE. NTOK .AND. ITOK .GT. 0) THEN
C
C---- Subroutine RDSCAL extracts from line, beginning at field ITOK,
C      sets of 3 arguments:- Label (= F1 or F2); Scale; Bfactor
C
C Input:  LSPRGI list of NLPRGI possible labels
C Output: ILPRGI  label number found
C         SCAL    scale
C         BB      temperature factor
C         ITOK    next item to read, = 0 if finished, .lt. 0 if error
C
I = ITOK
C                 ***************************************************
CALL RDSCAL(ITOK,LINE,IBEG,IEND,ITYPE,FVALUE,NTOK,NLPRGI,
+                  LSPRGI,ILPRGI,SCAL,BB)
C                 ***************************************************
IF (ITOK .EQ. I) THEN
C  Allow for old version of RDSCAL which doesn't reset ITOK
ITOK = ITOK + 3
ELSEIF (ITOK .EQ. -1) THEN
WRITE (LUNOUT, '(A)')
\$               ' **** No scale factor given ****'
KE = KE + 1
ELSEIF (ITOK .EQ. -2) THEN
WRITE (LUNOUT, '(A)')
\$               ' **** Unrecognized program label ****'
KE = KE + 1
ENDIF
C
IF (ITOK .GE. 0) THEN
SCL(ILPRGI) = SCAL
B(ILPRGI) = BB
COMPL(ILPRGI) = COMPLEX
WRITE (LUNOUT,FMT=6004),ILPRGI,
\$             LSPRGI(ILPRGI)(1:LENSTR(LSPRGI(ILPRGI))),
+             SCL(ILPRGI),B(ILPRGI),COMPL(ILPRGI)
6004           FORMAT (//
\$               '  Scale factor and Temperature factor for',7X,
\$               I4,A,4X,3F9.4)
GO TO 5
ENDIF
ENDIF
GO TO 20
C

ELSE IF(KEY.EQ.'END') THEN
WRITE(LUNOUT,'(/)')
GO TO 200
ELSE
Call ccperr(1,'  Key word unrecognised')
END IF
C
C----   End of keyword input
200  CONTINUE
C
IPRINT = 2
C     **********************************
CALL LROPEN(IFIL1,'HKLIN',IPRINT,IFAIL)
C     **********************************
IF(IFAIL.NE.0) THEN
CALL CCPERR(1, ' No input file???????')
END IF
C
C     ****************************************
CALL LRASSN(IFIL1,LSPRGI,NLPRGI,LOOKUP,CTPRGI)
C     *****************************************
C
C     =========================================
CALL LRCLAB(IFIL1,CLABS,CTUSRI,NCOL)
C     =========================================
C
C     Test correct column types.
NF1 = 0
NF2 = 0
DO 30 JDO30 = STARTFO,STARTFO+(NUMFO-1)*4,4
NF1 = NF1 + 1
NF2 = NF2 + LOOKUP(JDO30)
IF (LOOKUP(JDO30).NE.0 .NEQV. LOOKUP(JDO30+1).NE.0) THEN
WRITE (LUNOUT,FMT='(/A,I2,A,I2,A/)')
+    'Unmatched F',NF1,' or SIGF',NF1,' assigned.'
CALL CCPERR(1,' Unmatched F/SIGF pair')
ENDIF
IF (LOOKUP(JDO30+2).NE.0 .NEQV. LOOKUP(JDO30+3).NE.0) THEN
WRITE (LUNOUT,FMT='(/A,I2,A,I2,A/)')
+    'Unmatched PHIB',NF1,' or FOM',NF1,' assigned.'
CALL CCPERR(1,' Unmatched PHI/FOM pair')
ENDIF
IF (LOOKUP(JDO30+2).NE.0 .AND. LOOKUP(JDO30).EQ.0) THEN
WRITE (LUNOUT,FMT='(/A,I2,A,I2/)')
+    'PHI',NF1,' assigned but no corresponding F',NF1
CALL CCPERR(1,' Unmatched F/PHI pair')
ENDIF
30   CONTINUE

NF1 = 0
DO 40 JDO40 = STARTFC,STARTFC+(NUMFC-1)*2,2
NF1 = NF1 + 1
NF2 = NF2 + LOOKUP(JDO40)
IF (LOOKUP(JDO40).NE.0 .NEQV. LOOKUP(JDO40+1).NE.0) THEN
WRITE (LUNOUT,FMT='(/A,I2,A,I2,A/)')
+    'Unmatched FC',NF1,' or PHIC',NF1,' assigned.'
CALL CCPERR(1,' Unmatched FC/PHIC pair')
ENDIF
40    CONTINUE

NI1 = 0
NI2 = 0
DO 35 JDO35 = STARTI,STARTI+(NUMI-1)*2,2
NI1 = NI1 + 1
NI2 = NI2 + LOOKUP(JDO35)
IF (LOOKUP(JDO35).NE.0 .NEQV. LOOKUP(JDO35+1).NE.0) THEN
WRITE (LUNOUT,FMT='(/A,I2,A,I2,A/)')
+    'Unmatched I',NI1,' or SIGI',NI1,' assigned.'
CALL CCPERR(1,' Unmatched I/SIGI pair')
ENDIF
35    CONTINUE
C
C
IF ( NF2 .GT.0 .AND. NI2 .GT.0)
+      CALL CCPERR(1,' Cannot mix Fs and Is')

C
C     **********************
CALL LWOPEN(IFIL1,'HKLOUT')
C     **********************
C
C Append output label for FSUM and PHISUM
C     **********************
CALL LWASSN(IFIL1,LSPRGO,2,CTPRGO,1)
C     **********************
C
C---- Iflag = 0 means replace old title.
C
ICHK = 1
IFLAG = 0
IF (LTITLE) CALL LWTITL(IFIL1,TITLE,IFLAG)
C
C---- Loop for reflections----------------------------------------------------
C
10   CONTINUE
C
C     **************************
IF (ENDFIL) GO TO 100
CALL LRREFM(IFIL1,LOGMSS)
C     **************************
C
C       ICHK = MOD(NREF,1000)
FSUM = 0
PHISUM = 0.0
ASUM = 0
BSUM = 0.0
C     Test whether all terms there - if not will have to omit reflection
C
C
C     case (a): F, SIGF
DO 60 JDO60 = STARTFO,STARTFO+NUMFO*4-3,4
IF(LOOKUP(JDO60).EQ.0) GO TO 60
C     ...if SIGF already flagged as "Missing"
IF(LOGMSS(LOOKUP(JDO60)) .OR. LOGMSS(LOOKUP(JDO60+1)) .OR.
+ LOGMSS(LOOKUP(JDO60+2)) .OR. LOGMSS(LOOKUP(JDO60+3))) GO TO 10
IF(COMPL(JDO60).NE.0.0) PHI1 = PHI1 + COMPL(JDO60)*90
PHI1=DTOR*PHI1
ASUM = ASUM + F1*COS(PHI1)
BSUM = BSUM + F1*SIN(PHI1)
60   CONTINUE

C     case (c): FC, SIGFC
DO 70 JDO70 = STARTFC,STARTFC+NUMFC*2-2,2
IF(LOOKUP(JDO70).EQ.0) GO TO 70
IF(LOGMSS(LOOKUP(JDO70)) .OR. LOGMSS(LOOKUP(JDO70+1)))
+      GO TO 10
IF(COMPL(JDO70).NE.0.0) PHI1 = PHI1 + COMPL(JDO70)*90
PHI1=DTOR*PHI1
ASUM = ASUM + F1*COS(PHI1)
BSUM = BSUM + F1*SIN(PHI1)
70   CONTINUE

C     case (d): I, SIGI
DO 80 JDO80 = STARTI+1,STARTI+NUMI*2-1,2
IF(LOOKUP(JDO80).EQ.0) GO TO 80
IF(LOGMSS(LOOKUP(JDO80)) .OR. LOGMSS(LOOKUP(JDO80+1)))
+      GO TO 10
PHI1=0.0
ASUM = ASUM + F1*COS(PHI1)
BSUM = BSUM + F1*SIN(PHI1)
80   CONTINUE

FSUM = SQRT(ASUM*ASUM + BSUM*BSUM)
PHISUM=0.0
IF (FSUM.gt.0.0001) PHISUM=ATAN2(BSUM,ASUM)/DTOR
IF(ICHK.EQ.0)write(6,'(3i4,6f9.3)')in,ASUM,BSUM,FSUM,PHISUM
NREF = NREF + 1

C     ...Return for next reflection
GOTO 10

100  CONTINUE
CALL LRCLOS(IFIL1)
CALL LWCLOS(IFIL1,1)

WRITE(LUNOUT,'(//,A,/,A,//)')
+         ' CHANGES MADE TO MTZ FILE',
+         ' ========================'

WRITE(LUNOUT,'(/)')

CALL CCPERR(0,'Normal termination')
END

```