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 ============== PROGRAM RSTATS C ============== C C Original version written by S.E.V.Phillips C and subsequently modified by G.Fermi and A.C.Bloomer C C 16-Mar-1990 Peter Brick Imperial College C Program Converted to MTZ from LCF data file format and C many changes made to the organization of the code. C C Changes include:- C C 1) Temperature factor now applied before Rejection tests and C calculation of all R factors C C 2) Phase, S and Sigma Columns not required on input file C C 3) Introduce logarithmic scaling option C C 4) Remove option where a temperature factor C was refined but not applied C C 5) S recalculated if Different cell parameters specified C C 6) Calculation of Weights left in code - but Scale factor C removed from Weight calculation - I don't know where C these Weighting schemes come from! C C 7) Weighting Used in all Scale factor calculations C C 8) By default nothing Rejected and no reflections listed C C Still need to introduce better tests on convergence. C Some Numbers are calculated on the Fo Scale and others on the C output Scale. The r-factors can be Different if the B factor C is significant. DeltaF Rejection test - should this be done C with the data Scaled for output - or on the Fo Scale. C C---- Program description C =================== C C Scales two sets of F's together. C Calculates statistics for R on 1,2 and 3 Sigma data C and deletes bad data specified as Fc/Fo Ratio, C Sigma multiple, or /Fo-Fc/. C C .. Parameters .. INTEGER MaxResBin PARAMETER (MaxResBin=201) INTEGER MaxAmpBin PARAMETER (MaxAmpBin=21) INTEGER MaxZonBin PARAMETER (MaxZonBin=9) INTEGER MCOLS PARAMETER (MCOLS=500) C .. C .. Local Scalars .. REAL A11, A12, A22, Bneg, BShift, CorrelatCoeff, Dee, Dels, + DeltaFs, DETA, DF, DF2, DFDB, FC, FD, FO, LogFoFc, + OldTempFact, PHI, R, RB, RDEN, RecipSumCounter, + RecipSumWeights, RecipSumWgtFobs, Rfactor, Rfactor1Sigma, + Rfactor2Sigma, Rfactor3Sigma, RINC, RmsDeltaFsQ, + RmsWghtDelta, RNum, RT, S, S2, ScaleInitOLD, ScaleShift, + SFcalc, SIGFC, SIGFO, SMALL, SS, SSS, Sum1SigmaDelta, + Sum1SigmaFobs, Sum2SigmaDelta, Sum2SigmaFobs, Sum3SigmaDelta, + Sum3SigmaFobs, SumALLFobs, SumCounter, SumDeltaF, SumDiffs, + SumDiffsSquared, SumFcalc, SumFcalcSquared, SumFinalDiffs, + SumFobs, SumFobsFcalc, SumFobsSquared, SumWeights, + SumWghtDiffs, SumWghtFcalcSquared, SumWghtFobsFcalc, + SumWghtFobsSquared, TFAC, USED, V1, V2, W, WD2, WghtRfactor, + XSFO, XTEST, Y, SMAX INTEGER I, IFAIL, IFREE, IZON1, IZON2, IZON3, JCYC, JDO100, + JDO110, JDO20, JDO21, JDO30, JDO31,JDO80,JDO90, + KCurrAmpBin, KCurrResoBin,KCurrZonBin, KCycleCounter, + KFREE, LDUM, LUNIN, LUNOUT, MTZERROR, MTZPrintFlag, + NAmpBin, NAmpBin1, NLINES, NresBin, NresBin1, + NresoBinsUsed, NTEMP, NTotalRefls, Num1SigmaHKL, + Num2SigmaHKL, Num3SigmaHKL, NumDeltaReject, NumHKLINPUT, + NumHKLREAD, NumHKLSCALING, NumHKLUsed, NumMTZColSIN, + NumOUTPERCycle, NumRatioReject, NumResoReject, + NumSigmaReject, NumTotalReject, NumZEROF,NZonBin,NZonBin1 LOGICAL LASTCycleFlag, LCYCREJMONITOR, LNAN, LNAN1, MTZEOF, LSKIP, + LSTATSKIP CHARACTER DSTRNG*8, MTZVERSION*10, HSTRNG*80, IBKR*80 C .. C .. Local Arrays .. REAL AmpDeltaBins(MaxAmpBin,2), AmpFABSDiffBins(MaxAmpBin,2), + AmpAbsFcalcBins(MaxAmpBin,2), AmpAbsFobsBins(MaxAmpBin,2), + AmpFcalcBins(MaxAmpBin,2), AmpFCALSQBins(MaxAmpBin,2), + AmpFobsBins(MaxAmpBin,2), AmpFobsFCALBins(MaxAmpBin,2), + AmpFobsSQBins(MaxAmpBin,2), AmpWeightBins(MaxAmpBin,2), + AmpWghtDiffsBins(MaxAmpBin,2), + AmpWghtFobsSQUBins(MaxAmpBin,2), HKLIN(MCOLS), + HKLINSAVE(MCOLS), + RANGESMTZColS(2, MCOLS), ResoDeltaBins(MaxResBin,2), + ResoAbsFcalcBins(MaxResBin,2), ResoAbsFobsBins(MaxResBin,2), + ResoFABSDiffBins(MaxResBin,2), ResoFcalcBins(MaxResBin,2), + ResoFCALSQBins(MaxResBin,2), ResoFobsBins(MaxResBin,2), + ResoFobsFcalcBins(MaxResBin,2), + ResoFobsSQBins(MaxResBin,2), ResoWeightBins(MaxResBin,2), + ResoWghtDiffsBins(MaxResBin,2), + ResoWghtFobsSQUBins(MaxResBin,2), + ZonDeltaBins(MaxZonBin,2), ZonFABSDiffBins(MaxZonBin,2), + ZonAbsFcalcBins(MaxZonBin,2), ZonAbsFobsBins(MaxZonBin,2), + ZonFcalcBins(MaxZonBin,2), ZonFCALSQBins(MaxZonBin,2), + ZonFobsBins(MaxZonBin,2), ZonFobsFCALBins(MaxZonBin,2), + ZonFobsSQBins(MaxZonBin,2), ZonWeightBins(MaxZonBin,2), + ZonWghtDiffsBins(MaxZonBin,2), + ZonWghtFobsSQUBins(MaxZonBin,2) INTEGER IHKL(3), KAmpBinsCounter(MaxAmpBin,2), + KResoBinsCounter(MaxResBin,2), + KZonBinsCounter(MaxZonBin,2) LOGICAL LOGMSS(MCOLS) CHARACTER WORD(2)*8,HKLZON(8)*3,TITLE*70 C .. C .. External Functions .. INTEGER LENSTR EXTERNAL LENSTR C .. C .. External Subroutines .. EXTERNAL CCPDAT, CCPDPN, CCPERR, CCPFYP, CCPOPN, CCPRCS, + KEYIN, LRCLOS, LRINFO, LROPEN, LRREFL, LRREFM, + LRREWD, LWCLOS, LWHIST, LWREFL, MTZINI, WATE C .. C .. Intrinsic Functions .. INTRINSIC ABS, EXP, LOG, MAX, MIN, NINT, SQRT C .. C .. Common blocks .. COMMON /RSTTS/AbsFlag,KProcessFlag, KBKROUTFlag, LOUTPUTFlag, + LRatioRejectFlag, LSigmaRejectFlag, LDeltaRejectFlag, + LISTFlag, LWeightFlag, KPrintFlag, NFREE, NumFobsCol, + NumSIGFobsCol, NumFcalcCol, NumPHICol, NumSIGFcalCol, + NumFREECol, NumWeightCol, NColOUTSigma, NColOUTFcalc, + NColOUTSIGFCAL, NColOUTFREE, + NColOUTPHI, MTZFILENum, NumCycleS, RESMAX, RESMIN, RSBMAX, + RSBMIN, DMAX, DMIN, DSBMAX, DSBMIN, BinWidth, WidthInit, + WidthInit0, + ScaleInit, TempFactor, DiffPrintVALUE, DeltaRejectVALUE, + SigmaRejectVALUE, RejectRatio, WeightCoeffIC COMMON /RSTTSC/TITLE REAL BinWidth, DeltaRejectVALUE, DiffPrintVALUE, DMAX, DMIN, + DSBMAX, DSBMIN, RejectRatio, RESMAX, RESMIN, RSBMAX, RSBMIN, + ScaleInit, SigmaRejectVALUE, TempFactor, WidthInit, + WidthInit0 INTEGER KBKROUTFlag, KPrintFlag, KProcessFlag, LOUTPUTFlag, + LWeightFlag, MTZFILENum, NColOUTFcalc, NColOUTPHI, + NColOUTSIGFCAL, NColOUTFREE, NFREE, + NColOUTSigma, NumCycleS, NumFcalcCol, NumFobsCol, + NumFREECol, NumPHICol, NumSIGFcalCol, NumSIGFobsCol, + NumWeightCol LOGICAL LDeltaRejectFlag, LISTFlag, LRatioRejectFlag, + LSigmaRejectFlag,AbsFlag REAL WeightCoeffIC(4) DATA HKLZON/'EEE','EEO','EOE','EOO','OEE','OEO','OOE','OOO'/ C .. C .. Data statements .. DATA SMALL/1.0E-8/ C .. C NAmpBin = MaxAmpBin NresBin = MaxResBin NZonBin = MaxZonBin NAmpBin1 = NAmpBin - 1 NresBin1 = NresBin - 1 NZonBin1 = NZonBin - 1 LUNIN = 5 LUNOUT = 6 IFAIL = 0 LDUM = 80 SIGFO = 0.0 SIGFC = 0.0 PHI = 0.0 W = 0.0 LMON = 0 LMON1 = 0 LNAN = .TRUE. LNAN1 = .TRUE. C C---- LCycRejMonitor set true if too C many reflections Rejected on first Cycle C LCYCREJMONITOR = .FALSE. C C *********************************** CALL CCPFYP CALL MTZINI CALL CCPOPN(LUNIN, 'DATA', 5, 1, LDUM, IFAIL) CALL CCPOPN(LUNOUT, 'PRINTER', 6, 1, LDUM, IFAIL) CALL CCPRCS(LUNOUT, 'RSTATS', '$Date: 2002/08/06 14:13:45 $') CALL CCPDAT(DSTRNG) C ********************************** C C---- Print header banner C WRITE (LUNOUT, FMT=9000) 9000 FORMAT (/' RSTATS: Scale, Weights, AND R-FACTORS', + /' =====================================', //) C C---- Open input reflection file and Print brief header info C MTZFILENum = 1 MTZPrintFlag = 1 C C ************************************************ CALL LROPEN(MTZFILENum, 'HKLIN', MTZPrintFlag, MTZERROR) C ************************************************ C C---- Stop if error on opening input reflection file C IF (MTZERROR.NE.0) CALL CCPERR(1, ' ERROR IN LROPEN') C C---- Clear arrays for Amplitude RangesMtzCols and Resolution Bins C DO 20 JDO20 = 1, NresBin DO 10 JDO21 = 1, 2 KResoBinsCounter(JDO20, JDO21) = 0 ResoFABSDiffBins(JDO20, JDO21) = 0.0 ResoAbsFobsBins(JDO20, JDO21) = 0.0 ResoFobsBins(JDO20, JDO21) = 0.0 ResoWghtDiffsBins(JDO20, JDO21) = 0.0 ResoWghtFobsSQUBins(JDO20, JDO21) = 0.0 ResoDeltaBins(JDO20, JDO21) = 0.0 ResoFobsFcalcBins(JDO20, JDO21) = 0.0 ResoAbsFcalcBins(JDO20, JDO21) = 0.0 ResoFcalcBins(JDO20, JDO21) = 0.0 ResoFobsSQBins(JDO20, JDO21) = 0.0 ResoFCALSQBins(JDO20, JDO21) = 0.0 ResoWeightBins(JDO20, JDO21) = 0.0 10 CONTINUE 20 CONTINUE C DO 40 JDO30 = 1, NAmpBin DO 30 JDO31 = 1, 2 KAmpBinsCounter(JDO30, JDO31) = 0 AmpFABSDiffBins(JDO30, JDO31) = 0.0 AmpAbsFobsBins(JDO30, JDO31) = 0.0 AmpFobsBins(JDO30, JDO31) = 0.0 AmpWghtDiffsBins(JDO30, JDO31) = 0.0 AmpWghtFobsSQUBins(JDO30, JDO31) = 0.0 AmpDeltaBins(JDO30, JDO31) = 0.0 AmpFobsFCALBins(JDO30, JDO31) = 0.0 AmpAbsFcalcBins(JDO30, JDO31) = 0.0 AmpFcalcBins(JDO30, JDO31) = 0.0 AmpFobsSQBins(JDO30, JDO31) = 0.0 AmpFCALSQBins(JDO30, JDO31) = 0.0 AmpWeightBins(JDO30, JDO31) = 0.0 30 CONTINUE 40 CONTINUE C DO 41 JDO30 = 1, NZonBin DO 31 JDO31 = 1, 2 KZonBinsCounter(JDO30, JDO31) = 0 ZonFABSDiffBins(JDO30, JDO31) = 0.0 ZonAbsFobsBins(JDO30, JDO31) = 0.0 ZonFobsBins(JDO30, JDO31) = 0.0 ZonWghtDiffsBins(JDO30, JDO31) = 0.0 ZonWghtFobsSQUBins(JDO30, JDO31) = 0.0 ZonDeltaBins(JDO30, JDO31) = 0.0 ZonFobsFCALBins(JDO30, JDO31) = 0.0 ZonAbsFcalcBins(JDO30, JDO31) = 0.0 ZonFcalcBins(JDO30, JDO31) = 0.0 ZonFobsSQBins(JDO30, JDO31) = 0.0 ZonFCALSQBins(JDO30, JDO31) = 0.0 ZonWeightBins(JDO30, JDO31) = 0.0 31 CONTINUE 41 CONTINUE C C---- Read in key worded control info and set up C input and output reflection files C C ****** CALL KEYIN C ****** C C---- Use RANGE of FP to decide on sensible WidthInit. C IF (WidthInit0.EQ.1000.0) THEN C C ****************************************** CALL LRINFO(MTZFILENum, MTZVERSION, NumMTZColSIN, NumHKLINPUT, + RANGESMTZColS) C ****************************************** C WidthInit0 = RANGESMTZColS(2, NumFobsCol)/NAmpBin END IF WidthInit = WidthInit0 C WRITE (LUNOUT, FMT=9010) RSBMIN, RSBMAX, DSBMIN, DSBMAX, BinWidth 9010 FORMAT (/' Resolution limits for scaling: ', 2F10.3, ' Angstrom', + /' corresponding limits in 4(sin theta/lambda)**2: ', + 2F10.5, /' Bin Width: ', F9.4) WRITE (LUNOUT, FMT=9015) 9015 FORMAT (/' ALL reflections now written to output mtz file') WRITE (LUNOUT, FMT=9016) RESMIN, RESMAX 9016 FORMAT (/' Resolution limits for final statistics :', + 2F10.3,' (A)'/) C C---- Rejection criteria C WRITE (LUNOUT, FMT=9020) WidthInit0 9020 FORMAT (' Width of ranges on Fo before scaling: ', F10.1) IF (LISTFlag) WRITE (LUNOUT, FMT=9030) DiffPrintVALUE 9030 FORMAT (' Reflections listed if Scaled ABS(Fo-Fc).GT.', F7.1) WRITE (LUNOUT, FMT=9040) ScaleInit, TempFactor 9040 FORMAT (' Initial Scale factor (for Fc): K = ', F12.5, + /' Initial B factor: TempFactor = ', F10.5) IF (LSigmaRejectFlag) WRITE (LUNOUT, FMT=9050) SigmaRejectVALUE 9050 FORMAT (' Reflection Rejected if: Fo < SigmaRejectValue*Sig(Fo)', + /' where SigmaRejectValue = ', + F8.2) IF (LRatioRejectFlag) WRITE (LUNOUT, FMT=9060) RejectRatio 9060 FORMAT (' Reflection Rejected if: Fc/Fo < RejectRatio', + /' where RejectRatio = ', F8.2) IF (LDeltaRejectFlag) WRITE (LUNOUT, FMT=9070) DeltaRejectVALUE 9070 FORMAT (' Reflection Rejected if: Abs(Fo-Fc) > DeltaRejectValue' + , /' where DeltaRejectValue = ', + F7.0) WRITE (LUNOUT, FMT=9080) NumCycleS 9080 FORMAT (' Maximum Number Cycles of refinement: ', I5) C IF (KPrintFlag.EQ.0) THEN WRITE (LUNOUT, FMT=9090) 9090 FORMAT (' Print out only on Final Cycle') ELSE WRITE (LUNOUT, FMT=9100) 9100 FORMAT (' Print out on every Cycle') END IF C C---- Open file RSTATSBKR to output Final Scale and temperature factor C IF (KBKROUTFlag.EQ.1) CALL CCPDPN(13, 'RSTATSBKR', 'UNKNOWN', 'F', + 80, 0) C C---- Print Weighting info - need to expand this C C LWeightFlag C C 1 Sigma W=x1*(1/SD(FO))**2 C C 2 DELF W=1/(x1+x2*S) for S.GT.x4 (S=sintheta/lambda) C W=1/(x1+x2*S+x3*(x4-S)**2 for S.LT.x4 C C 3 DSIG "DELF" but multiplied by (1/SD(FO))**2 C C 4 EXP W=((1/SD(FO))**2)*x1/EXP(x2+x3*S) C IF (LWeightFlag.GT.0) WRITE (LUNOUT, FMT=9110) 9110 FORMAT (' Reflections Weighted by: ') IF (LWeightFlag.EQ.1) THEN WRITE (LUNOUT, FMT=9120) WeightCoeffIC(1) 9120 FORMAT (' W = A*(1/SigFo)**2 where A = ', F10.5) ELSE IF (LWeightFlag.EQ.2) THEN WRITE (LUNOUT, FMT=9130) WeightCoeffIC 9130 FORMAT (' W = 1/(A + B*s) for s.gt.C (s = sintheta/lambda)', + /' W = 1/(A + B*s + C(D-s)) for s.lt.D where', + /' A = ', F10.5, /' B = ', F10.5, + /' C = ', F10.5, /' D = ', F10.5) ELSE IF (LWeightFlag.EQ.3) THEN WRITE (LUNOUT, FMT=9140) WeightCoeffIC 9140 FORMAT (' W = 1/(A + B*s) * (1/SigFo)**2 for s.gt.C ', / + ' W = 1/(A + B*s + C(D-s)) * (1/SigFo)**2 fors.lt.D where' + , /' A = ', F10.5, /' B = ', F10.5, + /' C = ', F10.5, /' D = ', F10.5) ELSE IF (LWeightFlag.EQ.4) THEN WRITE (LUNOUT, FMT=9150) WeightCoeffIC(1), WeightCoeffIC(2), + WeightCoeffIC(3) 9150 FORMAT (' W = (1/SigFo)**2 * A/exp(B+C*s) where', + /' A = ', F10.5, /' B = ', F10.5, + /' C = ', F10.5) END IF C C---- Print output file info C IF (LOUTPUTFlag.EQ.0) THEN WRITE (LUNOUT, FMT=9160) 9160 FORMAT (' No output reflection file ') ELSE IF (LOUTPUTFlag.EQ.1) THEN WRITE (LUNOUT, FMT=9170) 9170 FORMAT (' Output file contains Fp,Fc,', + ' [SigFp,SigFc,phi,free,wt]') ELSE IF (LOUTPUTFlag.EQ.2) THEN WRITE (LUNOUT, FMT=9180) 9180 FORMAT (' Output file contains same Columns as input file') ELSE IF (KBKROUTFlag.EQ.1) THEN WRITE (LUNOUT, FMT=9190) 9190 FORMAT (/' B,K (for Fcalc), R-factor will be output (3F10.5)', + /' on to file RSTATSBKR') END IF C C---- Print Scale type info C IF (KProcessFlag.EQ.0) THEN WRITE (LUNOUT, FMT=9200) WRITE (LUNOUT, FMT=9210) 9200 FORMAT (' Scale and B determined by minimizing :', + /' Sum{ w(Fo-K*exp(-Bs)*Fc)**2 }') 9210 FORMAT (' Scale and B will be applied to Fcalc ', + /' (and to sig(Fcalc) if given)') ELSE IF (KProcessFlag.EQ.1) THEN WRITE (LUNOUT, FMT=9200) WRITE (LUNOUT, FMT=9220) 9220 FORMAT (' Scale and B will be applied to Fobs and sig(Fobs)') ELSE IF (KProcessFlag.EQ.2) THEN WRITE (LUNOUT, FMT=9200) WRITE (LUNOUT, FMT=9230) 9230 FORMAT (' Scale will be applied to Fobs, and B-factor to Fcalc' + ) ELSE IF (KProcessFlag.EQ.-1) THEN WRITE (LUNOUT, FMT=9240) 9240 FORMAT (' Scale determined by minimizing Sum{ (Fo-K**Fc)**2 }') WRITE (LUNOUT, FMT=9250) 9250 FORMAT (' Scale applied to Fobs and to SigFobs if given') ELSE IF (KProcessFlag.EQ.-2) THEN WRITE (LUNOUT, FMT=9240) WRITE (LUNOUT, FMT=9260) 9260 FORMAT (' Scale applied to Fcalc and to SigFcalc if given') ELSE IF (KProcessFlag.EQ.-3) THEN WRITE (LUNOUT, FMT=9270) 9270 FORMAT (' Scale and B determined by minimizing :', + /' Sum{ (log(Fo)-log(K*exp(-Bs)*Fc))**2 }') WRITE (LUNOUT, FMT=9210) ELSE IF (KProcessFlag.EQ.-4) THEN WRITE (LUNOUT, FMT=9270) WRITE (LUNOUT, FMT=9220) END IF C C---- End of set up C C Only one pass required if zero refinement Cycles requested C ie - just apply input Scale and temperature factors C IF (NumCycleS.EQ.0) THEN LASTCycleFlag = .TRUE. ELSE LASTCycleFlag = .FALSE. END IF C KCycleCounter = 1 C 50 CONTINUE C C---- Return here for each Cycle C ========================== C NumZEROF = 0 NumHKLREAD = 0 NumResoReject = 0 NumHKLSCALING = 0 NumOUTPERCycle = 0 NumSigmaReject = 0 NumRatioReject = 0 NumDeltaReject = 0 NresoBinsUsed = 0 C SumFobsFcalc = 0.0 SumFcalc = 0.0 SumAbsFcalc = 0.0 SumFobsSquared = 0.0 SumFcalcSquared = 0.0 SumDeltaF = 0.0 SumALLFobs = 0.0 SumALLAbsFobs = 0.0 SumFinalDiffs = 0.0 SumWeights = 0.0 SumDiffsSquared = 0.0 SumDiffs = 0.0 SumFobs = 0.0 SumAbsFobs = 0.0 SumCounter = 0.0 SumWghtDiffs = 0.0 SumWghtFobsSquared = 0.0 SumWghtFcalcSquared = 0.0 SumWghtFobsFcalc = 0.0 C Num1SigmaHKL = 0 Num2SigmaHKL = 0 Num3SigmaHKL = 0 C Sum1SigmaDelta = 0 Sum1SigmaFobs = 0 Sum1SigmaAbsFobs = 0 Sum2SigmaDelta = 0 Sum2SigmaFobs = 0 Sum2SigmaAbsFobs = 0 Sum3SigmaDelta = 0 Sum3SigmaFobs = 0 Sum3SigmaAbsFobs = 0 RT = 0 RB = 0 A11 = 0.0 A22 = 0.0 A12 = 0.0 V1 = 0.0 V2 = 0.0 RNum = 0.0 RDEN = 0.0 SMAX = 0.0 C C---- Rewind input file C IF (KCycleCounter.GT.1) CALL LRREWD(MTZFILENum) C IF (LASTCycleFlag) THEN C C---- For Final Cycle - turn on full Printing C KPrintFlag = 1 C C---- Print header for listed reflections C IF (LISTFlag) THEN WRITE (LUNOUT, FMT=9280) DiffPrintVALUE WRITE (LUNOUT, FMT=9290) 9280 FORMAT (' ++++ Reflections with /Fo - Fc/ .gt. ', F8.0, + ' ++++', / + ' h k l d(A) Fobs Fcalc Fo-Fc') 9290 FORMAT (' KEPT reflections are listed as Scaled for output') IF (LDeltaRejectFlag) WRITE (LUNOUT, FMT=9300) 9300 FORMAT (' RejectED reflections listed on Fo Scale') END IF END IF C C---- Start of reflection loop C CALL EQUAL_MAGIC(MTZFILENum,HKLIN,MCOLS) C 60 CONTINUE C LNAN = .TRUE. LNAN1 = .TRUE. C Set counter for whether reflection used or not.. JREF = 1 C Set flag for skipping scale calculation on last (ie output) cycle LSKIP = .FALSE. C Set flag for skipping statistics on last cycle LSTATSKIP = .FALSE. C---- Read reflection record. Items returned in file order C C *********************************** CALL LRREFL(MTZFILENum, SSS, HKLIN, MTZEOF) C *********************************** C IF (MTZEOF) GO TO 70 CALL LRREFM(MTZFILENum, LOGMSS) NumHKLREAD = NumHKLREAD + 1 C USED = -999 IHKL(1) = NINT(HKLIN(1)) IHKL(2) = NINT(HKLIN(2)) IHKL(3) = NINT(HKLIN(3)) C C---- Resolution s = 4*(sin(theta)/lambda)**2 = 1/d**2 C S = SSS S2 = 0.25*S C C---- Test reflection for Resolution limits C On last Cycle use RESO limits to reject from statistics C IF (LASTCycleFlag) THEN IF (SSS.LT.DMIN .OR. SSS.GT.DMAX) THEN LSTATSKIP = .TRUE. END IF ENDIF C IF (SSS.LT.DSBMIN .OR. SSS.GT.DSBMAX) THEN NumResoReject = NumResoReject + 1 IF (LASTCycleFlag) THEN LSKIP = .TRUE. ELSE GO TO 60 ENDIF END IF C IF (S.GT.SMAX) SMAX = S C C---- Assign input data to program variables. Test that required data C is present for this reflection. If some data missing then C increment NumZeroF Counter and skip to next reflection. C AsSume that first 3 Columns of reflection file contain indices C Output refl. on last cycle. C IF (LOGMSS(NumFobsCol) .OR. LOGMSS(NumFcalcCol)) THEN IF (.NOT. LSKIP) NumZEROF = NumZEROF + 1 IF (.NOT. LASTCycleFlag) GOTO 60 LSKIP = .TRUE. ELSE IF (ABS(HKLIN(NumFobsCol)).LT.0.0001 .OR. + ABS(HKLIN(NumFcalcCol)).LT.0.0001) THEN IF (.NOT. LSKIP) NumZEROF = NumZEROF + 1 IF (.NOT. LASTCycleFlag) GOTO 60 LSKIP = .TRUE. ENDIF ENDIF C IF (LOGMSS(NumFobsCol) ) THEN FO = 0.0 ELSE FO = HKLIN(NumFobsCol) IF(ABSFLAG) FO = ABS(FO) ENDIF IF (LOGMSS(NumFcalcCol) ) THEN FC = 0.0 ELSE FC = HKLIN(NumFcalcCol) IF(ABSFLAG) FC = ABS(FC) ENDIF C C---- Have to be careful that refl. are not counted twice for NumZEROF. C---- Assumed Fo=NaN then SIGFo=NaN - monitor first time a discrepancy found C IF (NumSIGFobsCol.GT.0) THEN IF (.NOT. LOGMSS(NumSIGFobsCol)) THEN IF (HKLIN(NumSIGFobsCol).LT.0.0001 .AND. ABS(FO).GT.0.0001) . THEN IF (.NOT. LSKIP) NumZEROF = NumZEROF + 1 IF (.NOT. LASTCycleFlag) GO TO 60 LSKIP = .TRUE. ENDIF SIGFO = HKLIN(NumSIGFobsCol) ELSE SIGFO = -10.0 IF (LOGMSS(NumSIGFobsCol) .NEQV. LOGMSS(NumFobsCol)) THEN LNAN = .FALSE. IF (LMON.EQ.0 ) THEN WRITE(6,'(/,A,A,/,A,A,/)') ' *** Warning: either ', . 'FP/SIGFP=MNF but other is not.',' This must mean the ', . 'column pair is inconsistent ***' LMON = 1 ENDIF ENDIF ENDIF END IF C IF (NumSIGFcalCol.GT.0) THEN IF (.NOT. LOGMSS(NumSIGFcalCol)) THEN IF (HKLIN(NumSIGFcalCol).LT.0.0001 .AND. ABS(FC).GT.0.0001) . THEN IF (.NOT. LSKIP) NumZEROF = NumZEROF + 1 IF (.NOT. LASTCycleFlag) GO TO 60 LSKIP = .TRUE. END IF SIGFC = HKLIN(NumSIGFcalCol) ELSE SIGFC = 0.0 IF (LOGMSS(NumSIGFcalCol) .NEQV. LOGMSS(NumFcalcCol)) THEN LNAN1 = .FALSE. IF (LMON1.EQ.0 ) THEN WRITE(6,'(/,A,A,/,A,A,/)') ' *** Warning: either ', . 'FC/SIGFC=MNF but other is not.',' This must mean the ', . 'column pair is inconsistent ***' LMON1 = 1 ENDIF ENDIF ENDIF END IF C IF(.NOT.LNAN .OR. .NOT. LNAN1) THEN IF (.NOT. LASTCycleFlag) GO TO 60 LSKIP = .TRUE. ENDIF C IFREE = 1 C IF (NumFREECol.GT.0) THEN IF (.NOT. LOGMSS(NumFREECol)) THEN C C---- Free Flag = 1 for included data; 2 for excluded. C IF (ABS(HKLIN(NumFREECol)-NFREE) .LE. 0.1) IFREE = 2 IF (IFREE.EQ.2) THEN IF (.NOT. LASTCycleFlag) GO TO 60 LSKIP = .TRUE. ENDIF ENDIF END IF C IF (NumPHICol.GT.0) PHI = HKLIN(NumPHICol) C*********** C C---- Apply Scale and temperature factor to Fc and form diffs - need to be careful! C TFAC = EXP(-TempFactor*S2) c USED = -999 IF(.NOT. LOGMSS(NumFcalcCol)) SFcalc = ScaleInit*FC*TFAC IF(AbsFlag) SFcalc = Abs(SFcalc) IF(.NOT. LOGMSS(NumFobsCol)) THEN FO = HKLIN(NumFobsCol) IF(AbsFlag) FO = Abs(FO) IF(.NOT. LOGMSS(NumFcalcCol)) THEN DeltaFs = FO - SFcalc RINC = ABS(DeltaFs) C Set this variable as a flag to save all this testing later! USED = 999 END IF END IF C C---- Skip all stats. if refl. not USED. Should only appear in NumZEROF refls. C IF (USED.GT.0.0) THEN C C---- R-factor as function of (Fo/SigFo) - before any Rejections C IF (LASTCycleFlag) THEN SumDeltaF = SumDeltaF + RINC SumALLFobs = SumALLFobs + FO SumALLAbsFobs = SumALLAbsFobs + ABS(FO) C IF (NumSigFobsCol .GT.0 .AND. .NOT. LOGMSS(NumSigFobsCol))THEN IF (ABS(FO).GE.SIGFO) THEN Num1SigmaHKL = Num1SigmaHKL + 1 Sum1SigmaDelta = Sum1SigmaDelta + RINC Sum1SigmaAbsFobs = Sum1SigmaAbsFobs + ABS(FO) Sum1SigmaFobs = Sum1SigmaFobs + FO C IF (ABS(FO).GE.2.0*SIGFO) THEN Num2SigmaHKL = Num2SigmaHKL + 1 Sum2SigmaDelta = Sum2SigmaDelta + RINC Sum2SigmaFobs = Sum2SigmaFobs + FO Sum2SigmaAbsFobs = Sum2SigmaAbsFobs + ABS(FO) C IF (ABS(FO).GE.3.0*SIGFO) THEN Num3SigmaHKL = Num3SigmaHKL + 1 Sum3SigmaDelta = Sum3SigmaDelta + RINC Sum3SigmaFobs = Sum3SigmaFobs + FO Sum3SigmaAbsFobs = Sum3SigmaAbsFobs + ABS(FO) END IF END IF END IF END IF ENDIF C C---- Test for weak reflection C IF (NumSigFobsCol .GT.0 .AND. .NOT.LOGMSS(NumSigFobsCol)) THEN IF (LSigmaRejectFlag .AND. + (ABS(FO).LT.SigmaRejectVALUE*SIGFO)) THEN NumSigmaReject = NumSigmaReject + 1 IF( .NOT. LASTCycleFlag) GO TO 60 LSKIP = .TRUE. END IF END IF C C---- Test for FC</) v Resln; ',a,':A:2,12:', + /': Fobs & Fcalc v Resln; ',a,':A:2,7,8: $$') WRITE (LUNOUT, FMT=9590) 9590 FORMAT (' Bin',/,' 1/resol^2',/, + ' Reso_in_Angstroms',/,' Refls_per_bin',/, + ' R_factor',/,' sqrt[(W*(Fo-Fc)^2)/(W*Fo^2)]',/, + ' Mean_AbsFobs',/,' Mean_AbsFc',/,' Mean_W*(Fo-Fc)^2',/, + ' (Fo-Fc)/(AbsFo)',/,' W_Mean_(Fo-Fc)',/, + ' ln_(Mean_AbsFo/Mean_AbsFc)',/, + ' Linear_Correlation_between_Fo_and_Fc $$',/, + ' Range 4(st/l)**2 d N R ', + ' WR / ', + 'Delsq ln(/) Correlation', /' $$') DO 80 JDO80 = 1, NresoBinsUsed IF (KResoBinsCounter(JDO80, JDO90).GT.0) THEN C C---- Bin Resolution is highest Resolution in Bin C S = JDO80*BinWidth SS = 1.0/SQRT(S) C C---- For highest resolution bin, set S to highest resolution reflection C IF (JDO80.EQ.NresoBinsUsed) THEN S = SMAX IF (S.GT.0.0) SS = 1.0/SQRT(S) ENDIF IF (ResoAbsFobsBins(JDO80, JDO90).EQ.0.0) THEN R = 0.0 WghtRfactor = 0.0 SumFinalDiffs = 0.0 CorrelatCoeff = 0.0 ELSE R = ResoFABSDiffBins(JDO80, JDO90)/ + ResoAbsFobsBins(JDO80, JDO90) SumFinalDiffs = ResoDeltaBins(JDO80, JDO90)/ + ResoAbsFobsBins(JDO80, JDO90) IF (ResoWghtFobsSQUBins(JDO80, JDO90).GT. + 0.0) WghtRfactor = SQRT(ResoWghtDiffsBins(JDO80, + JDO90)/ResoWghtFobsSQUBins(JDO80, JDO90)) SumCounter = KResoBinsCounter(JDO80, JDO90) FO = ResoAbsFobsBins(JDO80, JDO90)/SumCounter FC = ResoAbsFcalcBins(JDO80, JDO90)/SumCounter LogFoFc = LOG(FO/FC) RmsWghtDelta = SQRT(ResoWghtDiffsBins(JDO80, JDO90)/ + SumCounter) SumFobs = ResoFobsBins(JDO80, JDO90) SumFcalc = ResoFcalcBins(JDO80, JDO90) RmsDeltaFsQ = SQRT(ResoWghtDiffsBins(JDO80, JDO90)/ + ResoWeightBins(JDO80, JDO90)) XTEST = (ResoFobsSQBins(JDO80, JDO90)- + SumFobs*SumFobs/SumCounter)* + (ResoFCALSQBins(JDO80, JDO90)- + SumFcalc*SumFcalc/SumCounter) C IF (XTEST.LE.0.0) THEN CorrelatCoeff = -99.999 ELSE CorrelatCoeff = (ResoFobsFcalcBins(JDO80, JDO90)- + SumFobs*SumFcalc/SumCounter)/SQRT(XTEST) END IF END IF C WRITE (LUNOUT, FMT=9600) JDO80, S, SS, + KResoBinsCounter(JDO80, JDO90), R, WghtRfactor, FO, FC, + RmsWghtDelta, SumFinalDiffs, RmsDeltaFsQ, LogFoFc, + CorrelatCoeff 9600 FORMAT (1X, I5, F10.4, F10.2, I8, 9F10.3) C ResoFABSDiffBins(NresBin, JDO90) = ResoFABSDiffBins(NresBin, + JDO90) + ResoFABSDiffBins(JDO80, JDO90) ResoFobsBins(NresBin, JDO90) = ResoFobsBins(NresBin, + JDO90) + ResoFobsBins(JDO80, JDO90) ResoAbsFobsBins(NresBin, JDO90) = ResoAbsFobsBins(NresBin, + JDO90) + ResoAbsFobsBins(JDO80, JDO90) C KResoBinsCounter(NresBin, JDO90) = KResoBinsCounter(NresBin, C + JDO90) + KResoBinsCounter(JDO80, JDO90) ResoWghtDiffsBins(NresBin, JDO90) + = ResoWghtDiffsBins(NresBin, JDO90) + + ResoWghtDiffsBins(JDO80, JDO90) ResoWghtFobsSQUBins(NresBin, + JDO90) = ResoWghtFobsSQUBins(NresBin, JDO90) + + ResoWghtFobsSQUBins(JDO80, JDO90) ResoFcalcBins(NresBin, JDO90) = ResoFcalcBins(NresBin, + JDO90) + ResoFcalcBins(JDO80, JDO90) ResoAbsFcalcBins(NresBin, JDO90) = ResoAbsFcalcBins(NresBin, + JDO90) + ResoAbsFcalcBins(JDO80, JDO90) ResoDeltaBins(NresBin, JDO90) = ResoDeltaBins(NresBin, + JDO90) + ResoDeltaBins(JDO80, JDO90) ResoWeightBins(NresBin, JDO90) = ResoWeightBins(NresBin, + JDO90) + ResoWeightBins(JDO80, JDO90) ResoFobsFcalcBins(NresBin, JDO90) + = ResoFobsFcalcBins(NresBin, JDO90) + + ResoFobsFcalcBins(JDO80, JDO90) ResoFobsSQBins(NresBin, JDO90) = ResoFobsSQBins(NresBin, + JDO90) + ResoFobsSQBins(JDO80, JDO90) ResoFCALSQBins(NresBin, JDO90) = ResoFCALSQBins(NresBin, + JDO90) + ResoFCALSQBins(JDO80, JDO90) END IF 80 CONTINUE C IF (KResoBinsCounter(NresBin, JDO90).NE.0) THEN R = ResoFABSDiffBins(NresBin, JDO90)/ + ResoAbsFobsBins(NresBin, JDO90) SumCounter = KResoBinsCounter(NresBin, JDO90) WghtRfactor = SQRT(ResoWghtDiffsBins(NresBin, JDO90)/ + ResoWghtFobsSQUBins(NresBin, JDO90)) FO = ResoAbsFobsBins(NresBin, JDO90)/SumCounter FC = ResoAbsFcalcBins(NresBin, JDO90)/SumCounter LogFoFc = LOG(FO/FC) RmsWghtDelta = SQRT(ResoWghtDiffsBins(NresBin, JDO90)/ + SumCounter) SumFinalDiffs = ResoDeltaBins(NresBin, JDO90)/ + ResoFobsBins(NresBin, JDO90) SumFobs = ResoFobsBins(NresBin, JDO90) SumFcalc = ResoFcalcBins(NresBin, JDO90) RmsDeltaFsQ = SQRT(ResoWghtDiffsBins(NresBin, JDO90)/ + ResoWeightBins(NresBin, JDO90)) XTEST = (ResoFobsSQBins(NresBin, JDO90)- + SumFobs*SumFobs/SumCounter)* + (ResoFCALSQBins(NresBin, JDO90)- + SumFcalc*SumFcalc/SumCounter) CorrelatCoeff = (ResoFobsFcalcBins(NresBin, JDO90)- + SumFobs*SumFcalc/SumCounter)/SQRT(XTEST) C WRITE (LUNOUT, FMT=9610) 9610 FORMAT (' $$') WRITE (LUNOUT, FMT=9615) KResoBinsCounter(NresBin, JDO90), R, + WghtRfactor, FO, FC, RmsWghtDelta, SumFinalDiffs, + RmsDeltaFsQ, LogFoFc, CorrelatCoeff 9615 FORMAT (/' Overall Totals:', 10X, I8, 9F10.3, //) END IF 90 CONTINUE C C---- Amplitude ranges C DO 110 JDO110 = 1, KFREE IF (KAmpBinsCounter(NAmpBin, JDO110).EQ.0) GO TO 110 WRITE (LUNOUT, FMT=9620) 9620 FORMAT (/' Analysis against Fobs ', + /' (N.B. If the Fobs RANGE is narrow compared ', + /' with the expected Differences, the', + /' correlation Coefficient will be ', + /' essentially random)', /) WRITE (LUNOUT, FMT=9630) WORD(JDO110) 9630 FORMAT (' $TABLE :', a, ' R-Factor analysis against |Fobs| :') WRITE (LUNOUT, FMT=9640) WORD(JDO110), + TITLE(1:20),TITLE(1:20) 9640 FORMAT (' $GRAPHS: ', a, + ' R & Corrn v ; ',a,':N:5,3,10: ', + ':R v ; ',a,':N:6,3: $$') WRITE (LUNOUT, FMT=9650) 9650 FORMAT (' Bin',/,' Refls_per_bin',/,' R_factor',/, + ' sqrt[(W*(Fo-Fc)^2)/(W*Fo^2)]',/,' Mean_AbsFo',/, + ' Mean_AbsFc',/, + ' Mean_W*(Fo-Fc)^2',/,' (Fo-Fc)/(Fo)',/,' W_Mean_(Fo-Fc)',/, + ' Linear_Correlation_between_Fo_and_Fc $$',/, + ' Fo N R ', + ' WR / ', + ' Delsq Correlation', /' $$') DO 100 JDO100 = 1, NAmpBin1 IF (KAmpBinsCounter(JDO100, JDO110).GT.0) THEN IF (AmpAbsFobsBins(JDO100, JDO110).EQ.0.0) THEN R = 0.0 WghtRfactor = 0.0 SumFinalDiffs = 0.0 CorrelatCoeff = 0.0 ELSE R = AmpFABSDiffBins(JDO100, JDO110)/ + AmpAbsFobsBins(JDO100, JDO110) SumFinalDiffs = AmpDeltaBins(JDO100, JDO110)/ + AmpAbsFobsBins(JDO100, JDO110) IF (AmpWghtFobsSQUBins(JDO100, JDO110).GT. + 0.0) WghtRfactor = SQRT(AmpWghtDiffsBins(JDO100, + JDO110)/AmpWghtFobsSQUBins(JDO100, + JDO110)) SumCounter = KAmpBinsCounter(JDO100, JDO110) FO = AmpAbsFobsBins(JDO100, JDO110)/SumCounter FC = AmpAbsFcalcBins(JDO100, JDO110)/SumCounter RmsWghtDelta = SQRT(AmpWghtDiffsBins(JDO100, JDO110)/ + SumCounter) SumFobs = AmpFobsBins(JDO100, JDO110) SumFcalc = AmpFcalcBins(JDO100, JDO110) RmsDeltaFsQ = SQRT(AmpWghtDiffsBins(JDO100, JDO110)/ + AmpWeightBins(JDO100, JDO110)) XTEST = (AmpFobsSQBins(JDO100, JDO110)- + SumFobs*SumFobs/SumCounter)* + (AmpFCALSQBins(JDO100, JDO110)- + SumFcalc*SumFcalc/SumCounter) IF (XTEST.LE.0.0) THEN CorrelatCoeff = -99.999 ELSE CorrelatCoeff = (AmpFobsFCALBins(JDO100, JDO110)- + SumFobs*SumFcalc/SumCounter)/SQRT(XTEST) END IF END IF WRITE (LUNOUT, FMT=9660) JDO100, + KAmpBinsCounter(JDO100, JDO110), R, WghtRfactor, FO, FC, + RmsWghtDelta, SumFinalDiffs, RmsDeltaFsQ, CorrelatCoeff 9660 FORMAT (8X, I5, 5X, I8, 8F10.3) AmpFABSDiffBins(NAmpBin, JDO110) = AmpFABSDiffBins(NAmpBin, + JDO110) + AmpFABSDiffBins(JDO100, JDO110) AmpFobsBins(NAmpBin, JDO110) = AmpFobsBins(NAmpBin, + JDO110) + AmpFobsBins(JDO100, JDO110) AmpAbsFobsBins(NAmpBin, JDO110) = AmpAbsFobsBins(NAmpBin, + JDO110) + AmpAbsFobsBins(JDO100, JDO110) C KAmpBinsCounter(NAmpBin, JDO110) = KAmpBinsCounter(NAmpBin, C + JDO110) + KAmpBinsCounter(JDO100, JDO110) AmpWghtDiffsBins(NAmpBin, JDO110) + = AmpWghtDiffsBins(NAmpBin, JDO110) + + AmpWghtDiffsBins(JDO100, JDO110) AmpWghtFobsSQUBins(NAmpBin, JDO110) + = AmpWghtFobsSQUBins(NAmpBin, JDO110) + + AmpWghtFobsSQUBins(JDO100, JDO110) AmpFcalcBins(NAmpBin, JDO110) = AmpFcalcBins(NAmpBin, + JDO110) + AmpFcalcBins(JDO100, JDO110) AmpAbsFcalcBins(NAmpBin, JDO110) = AmpAbsFcalcBins(NAmpBin, + JDO110) + AmpAbsFcalcBins(JDO100, JDO110) AmpDeltaBins(NAmpBin, JDO110) = AmpDeltaBins(NAmpBin, + JDO110) + AmpDeltaBins(JDO100, JDO110) AmpWeightBins(NAmpBin, JDO110) = AmpWeightBins(NAmpBin, + JDO110) + AmpWeightBins(JDO100, JDO110) AmpFobsFCALBins(NAmpBin, JDO110) = AmpFobsFCALBins(NAmpBin, + JDO110) + AmpFobsFCALBins(JDO100, JDO110) AmpFobsSQBins(NAmpBin, JDO110) = AmpFobsSQBins(NAmpBin, + JDO110) + AmpFobsSQBins(JDO100, JDO110) AmpFCALSQBins(NAmpBin, JDO110) = AmpFCALSQBins(NAmpBin, + JDO110) + AmpFCALSQBins(JDO100, JDO110) END IF 100 CONTINUE C IF (KAmpBinsCounter(NAmpBin, JDO110).NE.0) THEN R = AmpFABSDiffBins(NAmpBin, JDO110)/ + AmpAbsFobsBins(NAmpBin, JDO110) SumCounter = KAmpBinsCounter(NAmpBin, JDO110) WghtRfactor = SQRT(AmpWghtDiffsBins(NAmpBin, JDO110)/ + AmpWghtFobsSQUBins(NAmpBin, JDO110)) FO = AmpAbsFobsBins(NAmpBin, JDO110)/SumCounter FC = AmpAbsFcalcBins(NAmpBin, JDO110)/SumCounter RmsWghtDelta = SQRT(AmpWghtDiffsBins(NAmpBin, JDO110)/ + SumCounter) SumFinalDiffs = AmpDeltaBins(NAmpBin, JDO110)/ + AmpFobsBins(NAmpBin, JDO110) SumFobs = AmpFobsBins(NAmpBin, JDO110) SumFcalc = AmpFcalcBins(NAmpBin, JDO110) RmsDeltaFsQ = SQRT(AmpWghtDiffsBins(NAmpBin, JDO110)/ + AmpWeightBins(NAmpBin, JDO110)) XTEST = (AmpFobsSQBins(NAmpBin, JDO110)- + SumFobs*SumFobs/SumCounter)* + (AmpFCALSQBins(NAmpBin, JDO110)- + SumFcalc*SumFcalc/SumCounter) CorrelatCoeff = (AmpFobsFCALBins(NAmpBin, JDO110)- + SumFobs*SumFcalc/SumCounter)/SQRT(XTEST) C C---- Overall Totals - Amplitude Bins C WRITE (LUNOUT, FMT=9670) 9670 FORMAT (' $$') WRITE (LUNOUT, FMT=9680) KAmpBinsCounter(NAmpBin, JDO110), R, + WghtRfactor, FO, FC, RmsWghtDelta, SumFinalDiffs, + RmsDeltaFsQ, CorrelatCoeff WRITE (LUNOUT, FMT=9690) 9680 FORMAT (/' Overall Totals:', 2X, I8, 8F10.3, //) 9690 FORMAT (/' --------------------------------------------', /) END IF 110 CONTINUE C C---- HKL zones C DO 130 JDO110 = 1, KFREE IF (KZonBinsCounter(NZonBin, JDO110).EQ.0) GO TO 130 WRITE (LUNOUT, FMT=9720) 9720 FORMAT (/' Analysis against HKL zones ', + /' (N.B. If there is pseudo symmetry', + ' it may distort this analysis.)',/) WRITE (LUNOUT, FMT=9730) WORD(JDO110) 9730 FORMAT (' $TABLE :', a, ' R-Factor analysis against |Zones| :') WRITE (LUNOUT, FMT=9740) WORD(JDO110), + TITLE(1:20) 9740 FORMAT (' $GRAPHS: ', a, + ' Corrn & R v Zone; ',a,':N:2,4,11: $$') WRITE (LUNOUT, FMT=9750) 9750 FORMAT (' Refl_Zones_Even_Odd',/, + ' Index',/,' Refls_per_bin',/,' R_factor',/, + ' sqrt[(W*(Fo-Fc)^2)/(W*Fo^2)]',/,' Mean_AbsFobs',/, + ' Mean_AbsFcalc',/,' Mean_W*(Fo-Fc)^2',/,' (Fo-Fc)/Fo',/, + ' W_Mean_(Fo-Fc)',/, + ' Linear_Correlation_between_Fo_and_Fc $$',/, + ' HKL_Zone IZ N R ', + ' WR / ', + ' Delsq Correlation', /' $$') DO 120 JDO100 = 1, NZonBin1 IF (KZonBinsCounter(JDO100, JDO110).GT.0) THEN IF (ZonAbsFobsBins(JDO100, JDO110).EQ.0.0) THEN R = 0.0 WghtRfactor = 0.0 SumFinalDiffs = 0.0 CorrelatCoeff = 0.0 ELSE R = ZonFABSDiffBins(JDO100, JDO110)/ + ZonAbsFobsBins(JDO100, JDO110) SumFinalDiffs = ZonDeltaBins(JDO100, JDO110)/ + ZonAbsFobsBins(JDO100, JDO110) IF (ZonWghtFobsSQUBins(JDO100, JDO110).GT. + 0.0) WghtRfactor = SQRT(ZonWghtDiffsBins(JDO100, + JDO110)/ZonWghtFobsSQUBins(JDO100, + JDO110)) SumCounter = KZonBinsCounter(JDO100, JDO110) FO = ZonAbsFobsBins(JDO100, JDO110)/SumCounter FC = ZonAbsFcalcBins(JDO100, JDO110)/SumCounter RmsWghtDelta = SQRT(ZonWghtDiffsBins(JDO100, JDO110)/ + SumCounter) SumFobs = ZonFobsBins(JDO100, JDO110) SumFcalc = ZonFcalcBins(JDO100, JDO110) RmsDeltaFsQ = SQRT(ZonWghtDiffsBins(JDO100, JDO110)/ + ZonWeightBins(JDO100, JDO110)) XTEST = (ZonFobsSQBins(JDO100, JDO110)- + SumFobs*SumFobs/SumCounter)* + (ZonFCALSQBins(JDO100, JDO110)- + SumFcalc*SumFcalc/SumCounter) IF (XTEST.LE.0.0) THEN CorrelatCoeff = -99.999 ELSE CorrelatCoeff = (ZonFobsFCALBins(JDO100, JDO110)- + SumFobs*SumFcalc/SumCounter)/SQRT(XTEST) END IF END IF WRITE (LUNOUT, FMT=9760) HKLZON(JDO100),JDO100, + KZonBinsCounter(JDO100, JDO110), R, WghtRfactor, FO, FC, + RmsWghtDelta, SumFinalDiffs, RmsDeltaFsQ, CorrelatCoeff 9760 FORMAT (8X,A3,3X, I5, 5X, I8, 8F10.3) ZonFABSDiffBins(NZonBin, JDO110) = ZonFABSDiffBins(NZonBin, + JDO110) + ZonFABSDiffBins(JDO100, JDO110) ZonFobsBins(NZonBin, JDO110) = ZonFobsBins(NZonBin, + JDO110) + ZonFobsBins(JDO100, JDO110) ZonAbsFobsBins(NZonBin, JDO110) = ZonAbsFobsBins(NZonBin, + JDO110) + ZonAbsFobsBins(JDO100, JDO110) C KZonBinsCounter(NZonBin, JDO110) = KZonBinsCounter(NZonBin, C + JDO110) + KZonBinsCounter(JDO100, JDO110) ZonWghtDiffsBins(NZonBin, JDO110) + = ZonWghtDiffsBins(NZonBin, JDO110) + + ZonWghtDiffsBins(JDO100, JDO110) ZonWghtFobsSQUBins(NZonBin, JDO110) + = ZonWghtFobsSQUBins(NZonBin, JDO110) + + ZonWghtFobsSQUBins(JDO100, JDO110) ZonFcalcBins(NZonBin, JDO110) = ZonFcalcBins(NZonBin, + JDO110) + ZonFcalcBins(JDO100, JDO110) ZonAbsFcalcBins(NZonBin, JDO110) = ZonAbsFcalcBins(NZonBin, + JDO110) + ZonAbsFcalcBins(JDO100, JDO110) ZonDeltaBins(NZonBin, JDO110) = ZonDeltaBins(NZonBin, + JDO110) + ZonDeltaBins(JDO100, JDO110) ZonWeightBins(NZonBin, JDO110) = ZonWeightBins(NZonBin, + JDO110) + ZonWeightBins(JDO100, JDO110) ZonFobsFCALBins(NZonBin, JDO110) = ZonFobsFCALBins(NZonBin, + JDO110) + ZonFobsFCALBins(JDO100, JDO110) ZonFobsSQBins(NZonBin, JDO110) = ZonFobsSQBins(NZonBin, + JDO110) + ZonFobsSQBins(JDO100, JDO110) ZonFCALSQBins(NZonBin, JDO110) = ZonFCALSQBins(NZonBin, + JDO110) + ZonFCALSQBins(JDO100, JDO110) END IF 120 CONTINUE C IF (KZonBinsCounter(NZonBin, JDO110).NE.0) THEN R = ZonFABSDiffBins(NZonBin, JDO110)/ + ZonAbsFobsBins(NZonBin, JDO110) SumCounter = KZonBinsCounter(NZonBin, JDO110) WghtRfactor = SQRT(ZonWghtDiffsBins(NZonBin, JDO110)/ + ZonWghtFobsSQUBins(NZonBin, JDO110)) FO = ZonAbsFobsBins(NZonBin, JDO110)/SumCounter FC = ZonAbsFcalcBins(NZonBin, JDO110)/SumCounter RmsWghtDelta = SQRT(ZonWghtDiffsBins(NZonBin, JDO110)/ + SumCounter) SumFinalDiffs = ZonDeltaBins(NZonBin, JDO110)/ + ZonFobsBins(NZonBin, JDO110) SumFobs = ZonFobsBins(NZonBin, JDO110) SumFcalc = ZonFcalcBins(NZonBin, JDO110) RmsDeltaFsQ = SQRT(ZonWghtDiffsBins(NZonBin, JDO110)/ + ZonWeightBins(NZonBin, JDO110)) XTEST = (ZonFobsSQBins(NZonBin, JDO110)- + SumFobs*SumFobs/SumCounter)* + (ZonFCALSQBins(NZonBin, JDO110)- + SumFcalc*SumFcalc/SumCounter) CorrelatCoeff = (ZonFobsFCALBins(NZonBin, JDO110)- + SumFobs*SumFcalc/SumCounter)/SQRT(XTEST) C C---- Overall Totals - HKL Zones C WRITE (LUNOUT, FMT=9770) 9770 FORMAT (' $$') WRITE (LUNOUT, FMT=9780) KZonBinsCounter(NZonBin, JDO110), R, + WghtRfactor, FO, FC, RmsWghtDelta, SumFinalDiffs, + RmsDeltaFsQ, CorrelatCoeff WRITE (LUNOUT, FMT=9790) 9780 FORMAT (/' Overall Totals:', 8X, I8, 8F10.3, //) 9790 FORMAT (/' --------------------------------------------', /) END IF 130 CONTINUE C IF (LOUTPUTFlag.GT.0) THEN C C---- Form one 80 char history line and write to output reflection file C It would be possible to include additional info on type C of scaling and to which Columns the scaling was applied C MTZPrintFlag = 1 NLINES = 1 WRITE (HSTRNG, FMT='(A,F9.3,A,F10.4,A,F10.4)') DSTRNG// + ' RSTATS R=', R, ' Scale =', ScaleInitOLD, ' Bfact = ', + OldTempFact C ******************************** CALL LWHIST(MTZFILENum, HSTRNG, NLINES) CALL LRCLOS(MTZFILENum) CALL LWCLOS(MTZFILENum, MTZPrintFlag) END IF IF (KBKROUTFlag .GT. 0) CLOSE(UNIT=13,STATUS='KEEP') C C ********************************************* CALL CCPERR(0, ' Normal Termination of rstats') C ********************************************* C END C C ================= SUBROUTINE KEYIN C ================= C C Routine sets default values and then reads in key worded control C information. C INTEGER MCOLS PARAMETER (MCOLS=500) INTEGER MSETS PARAMETER (MSETS=MCOLS) C .. Local Scalars .. INTEGER I, JAPPND, KWORK, LUNOUT, MTZAPPENDFlag, NLPRGI, NLPRGO, + NTOK, NumHKLINPUT, NumMTZColSIN LOGICAL LABOUTFlag, LEND CHARACTER KEY*4, MTZVERSION*10, TITLE*70, CWORK*400, + LABOUTSTORE*400, LINE*400 C .. C .. Local Arrays .. REAL FVALUE(200), RANGESMTZColS(2, MCOLS) INTEGER IBEG(200), IDEC(200), IEND(200), ITYP(200), + LOOKUP(MCOLS) CHARACTER CTPRGI(MCOLS)*1, CTPRGO(MCOLS)*1, + CVALUE(200)*4, LSPRGI(MCOLS)*30, + LSPRGO(MCOLS)*30 C C----Harvesting stuff C INTEGER NDATASETS,ISETS(MSETS),ISET,CSETID(MCOLS), + CSETOUT(MCOLS),SETID,IDUMMY CHARACTER*64 PNAME(MSETS),DNAME(MSETS), + PNAME_OUT(MCOLS),DNAME_OUT(MCOLS) C .. C .. External Functions .. INTEGER LENSTR EXTERNAL LENSTR C .. C .. External Subroutines .. EXTERNAL CCPERR, CCPUPC, GTPINT, GTPREA, LKYIN, LKYOUT, LRASSN, + LRINFO, LRRSOL, LWASSN, LWOPEN, LWTITL, PARSER, RDReso, + LRID, LRCLID, LWID, LWIDAS C .. C .. Intrinsic Functions .. INTRINSIC SQRT C .. C .. Common blocks .. COMMON /RSTTS/AbsFlag,KProcessFlag, KBKROUTFlag, LOUTPUTFlag, + LRatioRejectFlag, LSigmaRejectFlag, LDeltaRejectFlag, + LISTFlag, LWeightFlag, KPrintFlag, NFREE, NumFobsCol, + NumSIGFobsCol, NumFcalcCol, NumPHICol, NumSIGFcalCol, + NumFREECol, NumWeightCol, NColOUTSigma, NColOUTFcalc, + NColOUTSIGFCAL, NColOUTFREE, + NColOUTPHI, MTZFILENum, NumCycleS, RESMAX, RESMIN, RSBMAX, + RSBMIN, DMAX, DMIN, DSBMAX, DSBMIN, BinWidth, WidthInit, + WidthInit0, + ScaleInit, TempFactor, DiffPrintVALUE, DeltaRejectVALUE, + SigmaRejectVALUE, RejectRatio, WeightCoeffIC COMMON /RSTTSC/TITLE INTEGER KBKROUTFlag, KPrintFlag, KProcessFlag, LOUTPUTFlag, + LWeightFlag, MTZFILENum, NColOUTFcalc, NColOUTPHI, + NColOUTSIGFCAL, NColOUTFREE, NFREE, + NColOUTSigma, NumCycleS, NumFcalcCol, NumFobsCol, + NumFREECol, NumPHICol, NumSIGFcalCol, NumSIGFobsCol, + NumWeightCol LOGICAL LDeltaRejectFlag, LISTFlag, LRatioRejectFlag, + LSigmaRejectFlag,AbsFlag REAL WeightCoeffIC(4) C .. C .. Data statements .. DATA NLPRGI/9/ DATA LSPRGI/'H', 'K', 'L', 'FP', 'SIGFP', 'FC', 'PHIC', 'SIGFC', + 'FREE', 491*' '/ DATA CTPRGI/'H', 'H', 'H', 'F', 'Q', 'F', 'P', 'Q', 'I', + 491*' '/ DATA LOOKUP/-1, -1, -1, -1, 0, -1, 0, 0, 492*0/ DATA LSPRGO/'H', 'K', 'L', 'FP', 'SIGFP', 'FC', 'PHIC', 'SIGFC', + 'FREE', 'WT', 490*' '/ DATA CTPRGO/'H', 'H', 'H', 'F', 'Q', 'F', 'P', 'Q', 'I', + 'W', 490*' '/ DATA CSETID /MCOLS*0/ DATA CSETOUT/MCOLS*0/ C .. C LUNOUT = 6 TITLE = ' Output from RSTATS' NumCycleS = 6 RESMAX = -1.0 RSBMAX = RESMAX BinWidth = 0.010 WidthInit = 1000.0 LISTFlag = .FALSE. DiffPrintVALUE = 4000.0 WidthInit0 = 1000.0 ABSFLAG = .TRUE. LRatioRejectFlag = .FALSE. LSigmaRejectFlag = .FALSE. LDeltaRejectFlag = .FALSE. DeltaRejectVALUE = 99999.0 NFREE = 0 NumFobsCol = 0 NumFcalcCol = 0 NumSigFobsCol = 0 NumSigFcalCol = 0 NumPHICol = 0 NumFREECol = 0 C C---- Default for Process "FCAL" C KProcessFlag = 0 C C---- LabOutFlag set true when reflection file output C label editing required C LABOUTFlag = .FALSE. C C---- Default output NONE C KBKROUTFlag = 0 LOUTPUTFlag = 1 C C---- Default for Reject Ratio and Sigma C SigmaRejectVALUE = 0.0 RejectRatio = 0.0 C C---- Default for Scale and temperature factor C ScaleInit = 1.0 TempFactor = 0.0 C C---- Default for Print C KPrintFlag = 0 C C---- Defaults for WeightING_SCHEME NONE C LWeightFlag = 0 WeightCoeffIC(1) = 0.0 WeightCoeffIC(2) = 0.0 WeightCoeffIC(3) = 0.0 WeightCoeffIC(4) = 0.0 C C---- Start of loop to read control input C 10 CONTINUE KEY = ' ' NTOK = 200 LINE = ' ' C C *********************************************************** CALL PARSER(KEY, LINE, IBEG, IEND, ITYP, FVALUE, CVALUE, IDEC, + NTOK, LEND, .TRUE.) C *********************************************************** C IF (LEND) GO TO 50 C C---- Compare keyword with possible keys C IF (KEY.EQ.'TITL') THEN C C---- TITLE "string" C ===== C C The Key_Word is followed with the title information of up to C 80 characters C IF (NTOK.GE.2) TITLE = LINE(IBEG(2) :IEND(NTOK)) GO TO 10 c C---- Do not use absolute values of FP and FC C ============ C ELSE IF (KEY.EQ.'NOAB') THEN ABSFLAG = .FALSE. WRITE (LUNOUT,'(/,A,I8,/)') + ' Analysis done on Signed ( not absolute) values of FP and Fc' GO TO 10 C c C---- Free R value C ============ C ELSE IF (KEY.EQ.'FREE') THEN CALL GTPINT (2,NFREE,NTOK,ITYP,FVALUE) WRITE (LUNOUT,'(/,A,I8,/)') + ' Reflections treated as "free" if FreeRFlag = ',NFREE GO TO 10 C ELSE IF (KEY.EQ.'RESO') THEN C C---- ResoLUTION x1, x2 C ========== C C Read Dmin and Dmax in Angstrom (either order) C Default Dmin = 20.0 and Dmax = 2.0 C C If only one Number given then use this as Dmax C CALL RDReso(2, ITYP, FVALUE, NTOK, RESMAX, RESMIN, DMIN, DMAX) GO TO 10 ELSE IF (KEY.EQ.'RSCB') THEN C C---- RSCB x1, x2 C ==== C C Read Dmin and Dmax in Angstrom (either order) for scaling C Default Dmin = 20.0 and Dmax = 2.0 C C If only one Number given then use this as Dmax C CALL RDReso(2, ITYP, FVALUE, NTOK, RSBMAX, RSBMIN, DSBMIN, + DSBMAX) GO TO 10 ELSE IF (KEY.EQ.'SCAL') THEN C C---- Scale x C ===== C C ScaleInit Initial Scale factor C C ( Default is 1.0 ) C IF (NTOK.LE.1) THEN WRITE (LUNOUT, FMT=9000) ScaleInit 9000 FORMAT (//' ** No Argument given for Key_Word Scale', + /' Using default of ', F8.4, ' **') ELSE CALL GTPREA(2, ScaleInit, NTOK, ITYP, FVALUE) END IF GO TO 10 ELSE IF (KEY.EQ.'TEMP') THEN C C---- TEMPERATURE_FACTOR x C ================== C C TempFactor Initial temperature factor C C ( Default is 0.0 ) C IF (NTOK.LE.1) THEN WRITE (LUNOUT, FMT=9010) TempFactor 9010 FORMAT (// + ' ** No Argument given for Key_Word TEMPERATURE_FACTOR' + , /' Using default of ', F8.4, ' **') ELSE CALL GTPREA(2, TempFactor, NTOK, ITYP, FVALUE) END IF GO TO 10 ELSE IF (KEY.EQ.'WIDT') THEN C C---- Width_OF_Bins [one or both of RTHETA=x1, or FBinR=x2 ] C ============= C C RTHETA = x1 sets the Width of ranges of 4(sintheta/lambda)**2 C ( Default is 0.005 ) C C FBINR = x2 sets Width of ranges on F C ( Default is 1000.0 ) C IF (NTOK.LE.1) THEN WRITE (LUNOUT, FMT=9020) BinWidth, WidthInit0 9020 FORMAT (//' ** No arguments given for Key_Word Width_OF_Bins', + /' Using defaults of Rtheta = ', F8.4, ' and FBinr', + ' = ', F8.4, ' **') ELSE KWORK = 1 20 CONTINUE KWORK = KWORK + 1 IF (KWORK.LE.NTOK) THEN CWORK = LINE(IBEG(KWORK) :IEND(KWORK)) CALL CCPUPC(CWORK) IF (CWORK(1:4).EQ.'RTHE') THEN IF (NTOK.LE.KWORK) THEN WRITE (LUNOUT, FMT=9030) 9030 FORMAT (// + ' ** Key_Word input Error for < Width_OF_Bins RTHETA > **' + ) CALL CCPERR(1, ' ERROR IN RTHE') ELSE KWORK = KWORK + 1 CALL GTPREA(KWORK, BinWidth, NTOK, ITYP, FVALUE) GO TO 20 END IF ELSE IF (CWORK(1:4).EQ.'FBIN') THEN IF (NTOK.LE.KWORK) THEN WRITE (LUNOUT, FMT=9040) 9040 FORMAT (// + ' ** Key_Word input Error for < Width_OF_Bins FBINR > **' + ) CALL CCPERR(1, ' ERROR IN FBINR') ELSE KWORK = KWORK + 1 CALL GTPREA(KWORK, WidthInit0, NTOK, ITYP, FVALUE) GO TO 20 END IF END IF END IF END IF C GO TO 10 ELSE IF (KEY.EQ.'LIST') THEN C C---- LIST x (real value) [Optional] C ==== C Default is no listing! C ( Default is 4000.0 ) C C Sets value for listing of reflections with C /Fo - Scaled_Fc/ .gt. x C IF (NTOK.LE.1) THEN WRITE (LUNOUT, FMT=9050) DiffPrintVALUE 9050 FORMAT (//' ** No arguments given for Key_Word LIST', + /' Using Default of list hkl with Differences >', + F15.3, ' **') LISTFlag = .TRUE. GO TO 10 ELSE CALL GTPREA(2, DiffPrintVALUE, NTOK, ITYP, FVALUE) LISTFlag = .TRUE. END IF GO TO 10 ELSE IF (KEY.EQ.'CYCL') THEN C C---- CycleS x (integer value) [Optional] C ====== C C NumCycles Number of Cycles for scaling C C ( Default is 6 ) C IF (NTOK.LE.1) THEN WRITE (LUNOUT, FMT=9060) NumCycleS 9060 FORMAT (//' ** No arguments given for Key_Word CycleS', + /' Using Default Number of Cycles = ', I4, ' **') ELSE CALL GTPINT(2, NumCycleS, NTOK, ITYP, FVALUE) END IF GO TO 10 ELSE IF (KEY.EQ.'PRIN') THEN C C---- Print [one of either ALL or LAST ] C ===== C C ( Default is LAST ) C C ALL sets KPrintFlag on all Cycles C LAST sets KPrintFlag then Print out C on ONLY Final least squares Cycle C IF (NTOK.LE.1) THEN WRITE (LUNOUT, FMT=9070) 9070 FORMAT (//' ** No argument given for Key_Word Print - ', + ' Using default of LAST **') GO TO 10 ELSE CWORK = LINE(IBEG(2) :IEND(2)) CALL CCPUPC(CWORK) IF (CWORK(1:4).EQ.'LAST') THEN KPrintFlag = 0 ELSE IF (CWORK(1:3).EQ.'ALL') THEN KPrintFlag = 1 ELSE WRITE (LUNOUT, FMT=9080) CWORK(1:LENSTR(CWORK)) 9080 FORMAT (//' ** Warning for Key_Word Print,', + /' sub_key_word < ', A, ' > Not understood **') KPrintFlag = 0 END IF GO TO 10 END IF ELSE IF (KEY.EQ.'REJE') THEN C C---- Reject [one or more of Sigma=x1, Ratio=x2, Delta=x3 ] C ====== C C Ratio Sets value to Reject reflections C if K*TFAC*FC/FO .LT. RAT (ie. those with FC< **' + ) CALL CCPERR(1, ' ERROR IN SIGM') ELSE KWORK = KWORK + 1 CALL GTPREA(KWORK, SigmaRejectVALUE, NTOK, ITYP, FVALUE) LSigmaRejectFlag = .TRUE. GO TO 30 END IF ELSE IF (CWORK(1:4).EQ.'RATI') THEN IF (NTOK.LE.KWORK) THEN WRITE (LUNOUT, FMT=9110) 9110 FORMAT (// + ' ** Key_Word input Error for < Reject Ratio > **' + ) CALL CCPERR(1, ' ERROR IN RATI') ELSE KWORK = KWORK + 1 CALL GTPREA(KWORK, RejectRatio, NTOK, ITYP, FVALUE) LRatioRejectFlag = .TRUE. GO TO 30 END IF ELSE IF (CWORK(1:4).EQ.'DELT') THEN IF (NTOK.LE.KWORK) THEN WRITE (LUNOUT, FMT=9120) 9120 FORMAT (// + ' ** Key_Word input Error for < Reject Delta > **' + ) CALL CCPERR(1, ' ERROR IN DELT') ELSE KWORK = KWORK + 1 CALL GTPREA(KWORK, DeltaRejectVALUE, NTOK, ITYP, FVALUE) LDeltaRejectFlag = .TRUE. GO TO 30 END IF END IF END IF END IF C GO TO 10 ELSE IF (KEY.EQ.'OUTP') THEN C C---- OUTPUT [one of either C ====== C C NOHKL sets LOutputFlag = 0 C FOFC = 1 C C and optionally C C BKR sets KbkrOutFlag = 1 C C LOutputFlag = 0 no output file C = 1 output file has FP,SIGFP,FC,PHCAL [,Weight] C = 2 output file has same Columns as input [+ Weight] C C KbkrOutFlag = 1 Final B,K,R o/p on one line (3F10.5) on C file RSTATSBKR (ie, RSTATSBKR.DAT in default C directory unless otherwise assigned) C C (output file is Scaled as defined by Process) C IF (NTOK.LE.1) THEN WRITE (LUNOUT, FMT=9130) 9130 FORMAT (//' ** ERROR no argument given forKey_Word OUTPUT ', + '- Using default of FOFC **') ELSE KWORK = 1 40 CONTINUE KWORK = KWORK + 1 IF (KWORK.LE.NTOK) THEN CWORK = LINE(IBEG(KWORK) :IEND(KWORK)) C CALL CCPUPC(CWORK) C IF (CWORK(1:4).EQ.'NOHK') THEN IF(LOUTPUTFlag.NE.2) LOUTPUTFlag = 0 GO TO 40 ELSE IF (CWORK(1:4).EQ.'FOFC') THEN IF(LOUTPUTFlag.NE.2) LOUTPUTFlag = 1 GO TO 40 ELSE IF (CWORK(1:4).EQ.'ASIN') THEN WRITE (LUNOUT, FMT=9135) 9135 FORMAT(//,' ** ERROR: ASIN has been removed should use', + ' LABOUT ALLIN **') GO TO 40 ELSE IF (CWORK(1:3).EQ.'BKR') THEN KBKROUTFlag = 1 GO TO 40 ELSE WRITE (LUNOUT, FMT=9140) CWORK(1:LENSTR(CWORK)) 9140 FORMAT (//' ** Warning for Key_Word OUTPUT,', + /' sub_key_word < ', A, ' > Not understood **') GO TO 40 END IF END IF END IF GO TO 10 ELSE IF (KEY.EQ.'PROC') THEN C C---- Process [one of either C ======= C C FCAL sets KProcessFlag= 0 C FOBS = 1 C FOBC = 2 C SUMF = -1 C SUMC = -2 C LGFC = -3 C LGFO = -4 C C FCAL Apply Scale and B-factor to Fcalc and sigFc C Fobs Apply Scale and B-factor to Fo and sigFo C FOBC Apply Scale to Fo and sigFo, and B-factor to Fcalc C SUMF Calculate Scale by Sum(FoFc)/Sum(Fc**2) and C apply inverse of this to Fo C i.e. Temp factors are not refined and Scale C calculated without considering it C SUMC as "SumF" but apply Scale to Fc C LGFC Scale calculated by minimizing Sum((ln(Fo)-ln(Fc))**2) C LGFO Scale calculated same as LGFC but applied to Fo and SigFo C IF (NTOK.LE.1) THEN WRITE (LUNOUT, FMT=9150) 9150 FORMAT (//' ** No argument given for Key_Word Process - ', + ' Using default of Fcalc **') ELSE CWORK = LINE(IBEG(2) :IEND(2)) C C ************* CALL CCPUPC(CWORK) C ************* C IF (CWORK(1:4).EQ.'FCAL') THEN KProcessFlag = 0 ELSE IF (CWORK(1:4).EQ.'FOBS') THEN KProcessFlag = 1 ELSE IF (CWORK(1:4).EQ.'FOBC') THEN KProcessFlag = 2 ELSE IF (CWORK(1:4).EQ.'SUMF') THEN KProcessFlag = -1 ELSE IF (CWORK(1:4).EQ.'SUMC') THEN KProcessFlag = -2 ELSE IF (CWORK(1:4).EQ.'LGFC') THEN KProcessFlag = -3 ELSE IF (CWORK(1:4).EQ.'LGFO') THEN KProcessFlag = -4 ELSE WRITE (LUNOUT, FMT=9160) CWORK(1:LENSTR(CWORK)) 9160 FORMAT (//' ** Warning for Key_Word Process,', + /' sub_key_word < ', A, ' > Not understood **') KProcessFlag = 0 END IF GO TO 10 END IF ELSE IF (KEY.EQ.'WEIG') THEN C C---- WeightING_SCHEME [one of Sigma=x1, C ================ DELF=x1,x2,x3,x4 C DSIG=x1,x2,x3,x4 C EXP=x1,x2,x3 C or NONE ] C C ( Default is NONE ) C C Weight reflections according to one of the following schemes C C NONE No Weights to be applied C C Sigma W=x1*(K/SD(FO))**2 C C DELF W=1/(x1+x2*S) for S.GT.x4 (S=sintheta/lambda) C W=1/(x1+x2*S+x3*(x4-S)**2 for S.LT.x4 C C DSIG "DELF" but multiplied by (K/SD(FO))**2 C C EXP W=((K/SD(FO))**2)*x1/EXP(x2+x3*S) C IF (NTOK.LE.1) THEN WRITE (LUNOUT, FMT=9170) 9170 FORMAT (//, + ' ** No Argument given for Key_Word WeightING_SCHEME', + /,' Using default of NONE **') LWeightFlag = 0 ELSE KWORK = 2 CWORK = LINE(IBEG(KWORK) :IEND(KWORK)) CALL CCPUPC(CWORK) IF (CWORK(1:4).EQ.'NONE') THEN LWeightFlag = 0 GO TO 10 ELSE IF (CWORK(1:4).EQ.'SIGM') THEN IF (NTOK.LE.KWORK) THEN WeightCoeffIC(1) = 1.0 ELSE KWORK = KWORK + 1 CALL GTPREA(KWORK, WeightCoeffIC(1), NTOK, ITYP, FVALUE) END IF LWeightFlag = 1 GO TO 10 ELSE IF (CWORK(1:4).EQ.'DELF') THEN IF (NTOK.LE.KWORK) THEN WRITE (LUNOUT, FMT=9190) 9190 FORMAT (// + ' ** Key_Word input Error for < WeightING_SCHEME DELF > **' + ) CALL CCPERR(1, ' ERROR IN DELF') ELSE KWORK = KWORK + 1 CALL GTPREA(KWORK, WeightCoeffIC(1), NTOK, ITYP, FVALUE) KWORK = KWORK + 1 CALL GTPREA(KWORK, WeightCoeffIC(2), NTOK, ITYP, FVALUE) KWORK = KWORK + 1 CALL GTPREA(KWORK, WeightCoeffIC(3), NTOK, ITYP, FVALUE) KWORK = KWORK + 1 CALL GTPREA(KWORK, WeightCoeffIC(4), NTOK, ITYP, FVALUE) LWeightFlag = 2 GO TO 10 END IF ELSE IF (CWORK(1:4).EQ.'DSIG') THEN IF (NTOK.LE.KWORK) THEN WRITE (LUNOUT, FMT=9200) 9200 FORMAT (//, + ' ** Key_Word input Error for < WeightING_SCHEME DSIG > **' + ) CALL CCPERR(1, ' ERROR IN DSIG') ELSE KWORK = KWORK + 1 CALL GTPREA(KWORK, WeightCoeffIC(1), NTOK, ITYP, FVALUE) KWORK = KWORK + 1 CALL GTPREA(KWORK, WeightCoeffIC(2), NTOK, ITYP, FVALUE) KWORK = KWORK + 1 CALL GTPREA(KWORK, WeightCoeffIC(3), NTOK, ITYP, FVALUE) KWORK = KWORK + 1 CALL GTPREA(KWORK, WeightCoeffIC(4), NTOK, ITYP, FVALUE) LWeightFlag = 3 GO TO 10 END IF ELSE IF (CWORK(1:3).EQ.'EXP') THEN WeightCoeffIC(1) = 1.0 WeightCoeffIC(2) = 0.0 WeightCoeffIC(3) = 0.0 LWeightFlag = 4 KWORK = KWORK + 1 IF (NTOK.LT.KWORK) THEN WRITE(LUNOUT, FMT=9210) 9210 FORMAT(//,' Defaults employed for WEIGHT EXP') ELSE CALL GTPREA(KWORK, WeightCoeffIC(1), NTOK, ITYP, FVALUE) KWORK = KWORK + 1 IF (NTOK.GE.KWORK) + CALL GTPREA(KWORK, WeightCoeffIC(2), NTOK, ITYP, FVALUE) KWORK = KWORK + 1 IF (NTOK.GE.KWORK) + CALL GTPREA(KWORK, WeightCoeffIC(3), NTOK, ITYP, FVALUE) END IF GO TO 10 END IF END IF GO TO 10 ELSE IF (KEY.EQ.'LABI') THEN C C---- LABin Set up Column assignments for input reflection file C ===== C CALL LKYIN(1, LSPRGI, NLPRGI, NTOK, LINE, IBEG, IEND) GO TO 10 ELSE IF (KEY.EQ.'LABO') THEN C C---- LABOUT Set up Column assignments for output reflection file C ====== C C Store this line as we cannot edit labels until we know what C labels are required in the output C CWORK = LINE(IBEG(2):IEND(2)) CALL CCPUPC(CWORK) IF ( CWORK(1:5) .EQ. 'ALLIN') THEN LOUTPUTFlag = 2 DO 60 I = IBEG(2),IEND(2) LINE(I:I) = ' ' 60 CONTINUE ENDIF IF (LINE .NE. ' ') THEN LABOUTSTORE = LINE LABOUTFlag = .TRUE. ENDIF GO TO 10 C C---- END terminate input C === ELSE IF (KEY(1:3).EQ.'END') THEN GO TO 50 ELSE C C---- ???????????????? Unknown Key_Word C ================ C WRITE (LUNOUT, FMT=9220) KEY 9220 FORMAT (/' ** Key_Word Input of < ', A, ' >, Not Understood', + /' Ignoring this data line **') GO TO 10 END IF C 50 CONTINUE C C---- END Key_Word Input here C C Sort out Resolution limits C IF (RESMAX.LT.0.0) THEN C C---- No overall limits set, take them from file C CALL LRRSOL(MTZFILENum, DMIN, DMAX) RESMAX = 10000 RESMIN = 0.5 IF(DMIN.GT.0.0) RESMAX = 1.0/SQRT(DMIN) IF(DMAX.GT.0.0) RESMIN = 1.0/SQRT(DMAX) END IF C C---- No scaling limits set, set them = overall limits C IF (RSBMAX.LT.0.0) THEN RSBMAX = RESMAX RSBMIN = RESMIN DSBMIN = DMIN DSBMAX = DMAX END IF C C---- Set-up Column labels and Column types C for the input reflection file C C **************************************** CALL LRASSN(MTZFILENum, LSPRGI, NLPRGI, LOOKUP, CTPRGI) C **************************************** C C---- Dataset stuff C CALL LRID(MTZFILENum,PNAME,DNAME,ISETS,NDATASETS) C---- Get dataset IDs for columns and set up output IDs. IF (NDATASETS.GT.0) THEN CALL LRCLID(MTZFILENum,CSETID,IDUMMY) ENDIF C NumFobsCol = LOOKUP(4) NumSIGFobsCol = LOOKUP(5) NumFcalcCol = LOOKUP(6) NumPHICol = LOOKUP(7) NumSIGFcalCol = LOOKUP(8) NumFREECol = LOOKUP(9) C C---- Output reflection file C IF (LOUTPUTFlag.GT.0) THEN CALL LWOPEN(MTZFILENum, 'HKLOUT') C C---- Set-up Column labels and types depending on type of output C file requested C IF (LOUTPUTFlag.EQ.1) THEN C C---- Output h,k,l,s,Fp,SigFp,Fc,Sigfc,Phi,Weight,FreeR C We can only output Columns present in C input file + optional Weight C The array LookUp is Used to Flag those Columns C that must be present C MTZAPPENDFlag = 0 NumWeightCol = 0 I = 4 IF (NDATASETS.GT.0) + CSETOUT(I)=CSETID(LOOKUP(4)) C IF (NumSIGFobsCol.GT.0) THEN I = I + 1 NColOUTSigma = I LSPRGO(I) = 'SIGFP' CTPRGO(I) = 'Q' IF (NDATASETS.GT.0) + CSETOUT(I)=CSETID(LOOKUP(5)) END IF C IF (NumFcalcCol.GT.0) THEN I = I + 1 NColOUTFcalc = I LSPRGO(I) = 'FC' CTPRGO(I) = 'F' IF (NDATASETS.GT.0) + CSETOUT(I)=CSETID(LOOKUP(6)) END IF C IF (NumPHICol.GT.0) THEN I = I + 1 NColOUTPHI = I LSPRGO(I) = 'PHIC' CTPRGO(I) = 'P' IF (NDATASETS.GT.0) + CSETOUT(I)=CSETID(LOOKUP(7)) END IF C IF (NumSIGFcalCol.GT.0) THEN I = I + 1 NColOUTSIGFCAL = I LSPRGO(I) = 'SIGFC' CTPRGO(I) = 'Q' IF (NDATASETS.GT.0) + CSETOUT(I)=CSETID(LOOKUP(8)) END IF C IF (NumFREECol.GT.0) THEN I = I + 1 NColOUTFREE = I LSPRGO(I) = 'FREE' CTPRGO(I) = 'I' IF (NDATASETS.GT.0) + CSETOUT(I)=CSETID(LOOKUP(9)) END IF C IF (LWeightFlag.GT.0) THEN I = I + 1 NumWeightCol = I LSPRGO(I) = 'WT' CTPRGO(I) = 'W' C---- Weight gets same dataset ID as FP IF (NDATASETS.GT.0) + CSETOUT(I)=CSETID(LOOKUP(4)) END IF C NLPRGO = I C C---- Setup title on output file C JAPPND = 0 CALL LWTITL(MTZFILENum, TITLE, JAPPND) ELSE IF (LOUTPUTFlag.EQ.2) THEN C C---- Will write out all Columns in input file + [Weight] C MTZAPPENDFlag = 1 C IF (LWeightFlag.GT.0) THEN CALL LRINFO(MTZFILENum, MTZVERSION, NumMTZColSIN, + NumHKLINPUT, RANGESMTZColS) NumWeightCol = NumMTZColSIN + 1 NLPRGO = 1 LSPRGO(1) = 'WT' CTPRGO(1) = 'W' C---- Weight gets same dataset ID as FP IF (NDATASETS.GT.0) + CSETOUT(1)=CSETID(LOOKUP(4)) ELSE NumWeightCol = 0 NLPRGO = 0 END IF END IF C C---- Re-parse output label command that has been stored in LabOutStore C Note that if LOutputFlag=2 and C LWeightFlag=0 then C there are no program labels to edit C IF (LABOUTFlag) THEN CALL PARSER(KEY, LABOUTSTORE, IBEG, IEND, ITYP, FVALUE, + CVALUE, IDEC, NTOK, LEND, .TRUE.) CALL LKYOUT(1, LSPRGO, NLPRGO, NTOK, LABOUTSTORE, IBEG, IEND) END IF C C---- MtzAppendFlag =0 to replace labels, =1 to append labels and types C CALL LWASSN(MTZFILENum, LSPRGO, NLPRGO, CTPRGO, MTZAPPENDFlag) C Store the column dataset IDs: IF (NDATASETS.GT.0) THEN DO 115 JDO115 = 1,NLPRGO SETID = CSETOUT(JDO115) DO ISET = 1,NDATASETS IF (ISETS(ISET).EQ.SETID) THEN PNAME_OUT(JDO115) = PNAME(ISET) DNAME_OUT(JDO115) = DNAME(ISET) END IF ENDDO 115 CONTINUE CALL LWIDAS(MTZFILENum,NLPRGO,PNAME_OUT,DNAME_OUT, + MTZAPPENDFlag) ENDIF END IF C RETURN C END C C ===================================================== SUBROUTINE WATE(LWeightFlag, SIGFO, S, WeightCoeffIC, WK) C ===================================================== C C input parameters LWeightFlag C WeightCoeffic(1:4) which is read in as A,B,C,D C If LWeightFlag .GT. 0 then Fo and SigFo of all reflections C are Weighted according to one of the following schemes: C LWeightFlag=1 W=A*(K/SD(FO))**2 C LWeightFlag=2 W=1/(A+B*S) for S.GT.D (S=SINTHETA/LAMBDA) C W=1/(A+B*S+C*(D-S)**2 for S.LT.D C LWeightFlag=3 as 2 but multiplied by (K/SD(FO))**2 C LWeightFlag=4 W=((K/SD(FO))**2)*A/EXP(B+C*S) C C .. Scalar Arguments .. REAL S, SIGFO, WK INTEGER LWeightFlag C .. C .. Array Arguments .. REAL WeightCoeffIC(4) C .. C .. Intrinsic Functions .. INTRINSIC EXP C .. C IF (LWeightFlag.EQ.1) THEN C C---- Absolute Weights 1/(Sigma f )**2 C WK = WeightCoeffIC(1)/ (SIGFO*SIGFO) ELSE IF (LWeightFlag.EQ.2 .OR. LWeightFlag.EQ.3) THEN IF (S.GE.WeightCoeffIC(4)) THEN WK = 1.0/ (WeightCoeffIC(2)*S+WeightCoeffIC(1)+ + ((WeightCoeffIC(4)-S)**2)*WeightCoeffIC(3)) IF (LWeightFlag.NE.2) WK = WK/ (SIGFO*SIGFO) ELSE WK = 1.0/ (WeightCoeffIC(2)*S+WeightCoeffIC(1)) IF (LWeightFlag.NE.2) WK = WK/ (SIGFO*SIGFO) END IF ELSE IF (LWeightFlag.EQ.4) THEN WK = WeightCoeffIC(3)*S + WeightCoeffIC(2) WK = WeightCoeffIC(1)/ (SIGFO*SIGFO*EXP(WK)) END IF C RETURN C END