C C************************************************************* C C Program Converted to MTZ from LCF data file format C C Date: Tue Nov 19 12:24:08 GMT 1991 C by: Eleanor Dodson (York University) C C The development of the MTZ files and associated software mark 1 C is part of the masterplan of the ESF/EACBM Working Group 2.1 for C better Protein Crystallographic software for Europe. C C************************************************************** C C C =============== PROGRAM SIGMAAM C =============== C C C Changes made to SIGMAA.FOR C - Kevin Cowtan, 18/6/91 C C When combining with MIR data the Hendrickson and C Lattmann coefficients for the MIR data are normally C supplied in the input intensity file. If they are not C given however then they may be calculated from the MIR phase C and figure of merit by the relationships C C X = SIM^-1 ( FOM ) (acentric) C X = ATANH ( FOM ) (centric) C A = X COS ( PHI ) C B = X SIN ( PHI ) C C = 0 C D = 0 C where SIM(X) = I1(X) / I0(X) C Note, 0 < FOM < 1, 0 < X < infinity C C In the current version of the program a lookup table is generated C (ABTAB) which contains SIM(X) for 20 values of X. This is then used C as a direct lookup for the values of X from the FOM. Thus the C calculation X = SIM ( FOM ) is performed, which is the inverse of C the desired function. In the corrected version a table lookup and C interpolation is performed to perform the inverse function C X = SIM^-1 ( FOM ). The correct function for centric reflection C is also included. C C The problem showed up in the output of the program in that the mean MIR C figures of merit listed are not consistent with those in the input data. C This has been cured, and the new phases and figures of merit have been C shown to be better than before. C C Note that in the case where H-L coefficients are provided in the input C file the problem does not arise and there should be no change in the C results. C C C SIGMAA -- IMPROVED FOURIER COEFFICIENTS USING CALCULATED PHASES C R. J. READ C C----- HISTORY ------------------------------------------------ C C March 1995 Eleanor Dodson, Adam Ralph: C Determined effort to get back into step with Randy C Now outputs Randys coefficients for "2mFo-DFc" type maps C and "mFo-DFc" type maps. C C Option to read in two sets of HL cooefficents and ADD C C Two options to modify SIGMAA values. C Randy prefers to read in set obtained from somewhere else: C Victor Lamsin wants a "damping" factor. C C January 1991 Eleanor Dodson: replace FOM with W . C Add correlations to Phase ananysis. C September 1989 Phil Evans MRC LMB C Some tidying up, change of keywords, different output C columns. Options 1 & 2 combined to give both difference C & 2Fo-Fc type coefficients at the same time C C September 1988: E.J. Dodson: When using IOPT=1 which should give C terms for a weighted difference fourier, C the program assumes that FO and FC are C on the same scale; this is not necessary C and is not assumed when IOPT=2 OR 3. C C Now it calculates C C (m*EO - Ec) and converts these terms back to the FO scale. C C ie SQRT(EPS*SIGMAN)*(M*EOBS-ECALC),ALPHA CALC, C C February 1988: E.J. Dodson : (1) Modifications to read intensity files. C discarded dumper C (2) Input Through call cards C (3) Centric and epsilon zones found C by using symmetry for appropriate C spacegroup. C (4) Structure factors read from one C intensity file containing Fobs and C any FC from partial structures C which may be needed. C (5) Output to intensity file of new weights C M and D and if requested C a combined phase. C (6) The LWREFL call writes out the C contents of the input file, C plus these items:- C (7) cell and star only used to C calculate sin limits - C not needed for intensity C move dlim(..) to /tabs/ C COMMON /CELPAR/ cell(6),star(6),dlim(51) C C C (8) NSTEP is number of steps for C phase probability integration. C must be set to .le. dimension C of COSA and SINA. C C August 1986: R.J. Read: LAST MODIFICATIONS C C C Randy's notes - plus extensions C C There are three(now 4) options chosen through IOPT(ejd). C Coefficients are ouput for different types of maps: C EJD IOPT values are: IOPT = 1, partial FCs , no MIR info C IOPT = 3, MIR info, plus partial FCs or 2nd MIR set C Look for NPS MIR2 C IOPT = 4, Do PHANAL only C C---- 1a. FO-FC type map coeffs -- these are of the form C (M*FOBS - D*FCALC, ALPHA CALC) for partial phases C 1a is Randy's IOPT 1 C C Combined phase difference map coeffs. C C Centric and acentric data: C 1b. (Mcomb*FOBS, ALPHA comb - D*FCALC(1), ALPHA CALC(1) C (M,D is described under 2) C 1b is Randys IOPT 4 C C---- 2a. Reduced-bias Native map coefficients. C (similar to 2Fo-Fc type map coeffs)-- These are of the form C Centric data: C M*FOBS, ALPHA CALC, C Acentric data: C 2*M*FOBS-D*FCALC, ALPHA CALC, C where D is as defined by Luzzati (1952) and C SIGMAA is as defined by R. Sinivasan (Acta Cryst. 20: 143(1966)). C EPS is the epsilon value, ie. multiplicity for this zone C (see Rollett, p. 126). C SIGMAN = , SIGMAP = C C M = = FIGURE OF MERIT, C which is calculated from EOBS, ECALC and SIGMAA. C one should note that SIGMAA = D*SQRT(SIGMAP/SIGMAN) C and that, eg., EOBS = FOBS/SQRT(EPS*SIGMAN). C the map coeffs arise from the approximation that C M*EOBS,ALPHA CALC = 0.5*EOBS,ALPHA TRUE + C 0.5*SIGMAA*ECALC,ALPHA CALC C (R. J. READ, ACTA CRYST. A42: 140 (1986)) C 2a is Randys IOPT 2 C C---- 2b. Combined phase map coeffs modified to reduce bias. C Centric data: C Mcomb*FOBS, ALPHA comb, C Acentric data: C For mir phase info, M*|FO|, ALPHA MIR = FO + DELTA C (where delta is a random error). C For partial structure, C M*|FO|, ALPHA PAR = FO/2 + D*FC/2 + DELTA C Assume that combined phase will give a value that varies C between these as a linear function of w, which varies from C 0 to 1 with amount of partial structure influence. C We can include several partial structures to get: C M*|FO|, ALPHA COMB = ((2-SUM(Wi))/2)*FO + SUM((Wi*Di/2)*FCi) C (ignoring by now the random error component) C Then solve for fo (remember these are vectors with phases): C FO = (M*|FO|,ALPHA COM - SUM((Wi*Di/2)*FCi)) / (1-SUM(Wi)/2) C For the weight, we use the variation of information for the C phase probability from that source, divided by the sum of C the VARINF's from all of the sources (so the W's add up to 1.) C 2b is Randys IOPT 3 C c Negative iopt. C This is NOT an option in the MTZ version EJD C C Kept here for information only.... C C If IOPT is given as a negative number, all the relevant C numbers from all the sources of phase information are DUMPed C on unit 8, so that they can be used to phase delta F's C or to calculate any desired function of these numbers. C the header for the DUMP file is created by s/r header C (see there for explanation). for IOPT of 1 or 2 (calculated C phases only), reflection info is DUMPed by s/r sfout. C for iopt of 3 (combined phases), s/r DUMPer is used. C (look there for the format; the array rDUMP is explained in C s/r combin) C C rDUMP contains, eventually: C 1) combined phase (mir and partial structure, if available), C 2) combined figure of merit, 3) combined variation of C information (varinf); then for each partial structure, C 1) FC, 2) PHASE, 3) W, 4) VARINF; ie. up to 12 slots for C partial structure stuff; then for mir data, C 1-4) hendrickson-lattman coefficients a, b, c, d, C 5) phase, 6) W, 7) varinf C varinf slots of rDUMP set to -1.0 so that, if not C filled in, this will indicate no info from this source C C from s/r mirfil C ioff=4*nps C Rdump(IOFF+1) = A C Rdump(IOFF+2) = B C Rdump(IOFF+3) = C C Rdump(IOFF+4) = D C Rdump(IOFF+5) = ALF C Rdump(IOFF+6) = W C C FOR VAX/VMS, EXPLICITLY OPEN UNIT 8 TO ALLOW LONG LINES. C C IF (dump) OPEN(UNIT=8,STATUS='NEW',RECL=255) C C********************* END OF HISTORY *************************** C C INTEGER ISIZ PARAMETER (ISIZ=1000000) C .. Scalars in Common .. REAL ARGMAX,CONV,DELST2,SMAX,ST2MIN,WMAX,WMAX2 INTEGER IERR,IOPT,ISYSR,ISYSW,NBIN,NCENT,NCOL,NEZONE,NMON,NPS, + NSTEP,MIR2,IRFSIG C .. C .. Arrays in Common .. REAL COSA,CPROJ,DLIM,EPZONE,ETAB,SIGMAA,SIGMAN,SIGMAP,SINA,ACHANG C .. C .. Local Scalars .. INTEGER IPS C .. C .. External Subroutines .. EXTERNAL CCPERR,CCPFYP,COMBIN,DETSIG,LRCLOS,LRREWD,LWCLOS,MTZINI, + RCARDS,SFOUT,TABLES C .. C .. Common blocks .. COMMON /CNTROL/NBIN,ST2MIN,DELST2,NPS,MIR2,IOPT,SMAX,IERR,WMAX, + WMAX2,SMIN,SMINCB,SMAXCB,IRFSIG COMMON /CP/CPROJ(3,20),NCENT COMMON /EPS/EPZONE(4,20),NEZONE COMMON /IO/ISYSR,ISYSW,NCOL,CONV,NMON COMMON /SIGS/SIGMAA(3,50),SIGMAN(3,50),SIGMAP(3,50),ACHANG(3) COMMON /TABS/NSTEP,COSA(120),SINA(120),ETAB(4001),ARGMAX,DLIM(51) C .. C C ****** CALL CCPFYP call ccp4h_init() call ccp4h_summary_beg() call ccp4h_header('SIGMAA','sigmaa',1) call ccp4h_summary_end() call ccp4h_pre_beg() CALL CCPRCS(6,'SIGMAA','$Date: 2002/09/27 09:53:07 $') CALL MTZINI C ****** C NSTEP = 120 CONV = 3.14159/180.0 C C---- Read input and control cards. C C ****** CALL RCARDS C ****** C C----Ejd 1-sep-1988 We need to call detsig to get sigman(1,*) C even if nps=0 C sigman(...) is sum(fp**2) v sin theta.. C IPS = 1 IF (NPS.EQ.0) IPS = 0 10 CONTINUE C IF (IOPT.EQ.3) WRITE (ISYSW,FMT=6000) IPS 6000 FORMAT (' SIGMAA Determination for Partial Structure',I2) C C---- Determine sigmaa for each range of sin(theta)/lambda of the C calculated structure factors for use in figure of merit C calculation. C C *********** CALL DETSIG(IPS) C *********** C C---- MTZ rewind C C ********* CALL LRREWD(1) C ********* C IPS = IPS + 1 IF (IPS.LE.NPS) GO TO 10 C C---- If phases are not to be combined, write out appropriate map C coefficients for the calculated structure. C IF (IOPT.LT.3) THEN C C ***** CALL SFOUT C ***** C ELSE C C---- For option 3, combine hendrickson-lattman coeffs and C calculate phases, figures of merit and information content by C numerical integration. C C---- Set up cosine, sine and exp tables C C ******* CALL TABLES C ******* C C---- Calculate hendrickson-lattman coefficients for C partial structure(s), then produce fourier coefficients C using mir phase information (in form of h-l coeffs) and C partial structures. C C ****** CALL COMBIN C ****** C END IF C C---- Close input and output MTZ files. C C ********************************************** CALL LRCLOS(1) CALL LWCLOS(1,1) CALL CCPERR(0, 'Normal termination') C ********************************************** C END C C ================= SUBROUTINE RCARDS C ================= C C .. Parameters .. INTEGER MCOLS PARAMETER (MCOLS=500) INTEGER NKEYS,NPAR PARAMETER (NKEYS=11,NPAR=200) C .. C .. Scalars in Common .. REAL ARGMAX,CONV,DELST2,SMAX,ST2MIN,WMAX,WMAX2,SMIN,SMINCB,SMAXCB INTEGER IAC1,IAC2,IAC3,IAO,IAO2,ICOA,ICOB,ICOC,ICOD,ICHK,IERR, + IFC1, + IFC2,IFC3,IFO,IFS,IOPT,IPRINT,ISYSR,ISYSW,IW,IW2,JX,JY,JZ, + LASYM,LATOMS,LPRINT,MATOMS,NBIN,NCOL,NMON, + NPS,NSG,NSTEP,NSYM,MIR2,IRFSIG C .. C .. Arrays in Common .. REAL COSA,DLIM,ETAB,PERM,RSYM,RSYMT,SIGMAA, + SIGMAN,SIGMAP,SINA C .. C .. Local Scalars .. REAL RMAX,RMIN,SLIM,STMAX,STMIN,STOL2,RMINCB,RMAXCB INTEGER I,IFAIL,ITOK,LDUM, NBIN1,NLPRGI,NLPRGO,NMULT,NSPGRP, + NSYMP,NTOK,NSMP,N,ICHANG,MPRINT LOGICAL LEND CHARACTER KEY*4,NAMSPG*10,PGNAME*10,TITLE*80,LABOUT_SAVE*600, + LINE*600,LTYPE*1 C .. C .. Local Arrays .. REAL CELL(6),FVALUE(NPAR) INTEGER IBEG(NPAR),IDEC(NPAR),IEND(NPAR),ITYP(NPAR),LOOKUP(MCOLS) CHARACTER CTPRGI(MCOLS)*1,CTPRGO(MCOLS)*1,CTYPS(MCOLS)*1, + CVALUE(NPAR)*4,KEYWRD(NKEYS)*4,CLABS(MCOLS)*30, + LSPRGI(MCOLS)*30,LSPRGO(MCOLS)*30 C .. C .. External Subroutines .. EXTERNAL CCPDPN,CCPERR,CCPUPC,CENTRIC,EPSLN,GTPINT,LKYIN, + LKYOUT,LRASSN,LRCELL,LRCLAB,LROPEN,LWASSN,LWOPEN,LWTITL, + PARSER,RDRESO,RDSYMM, LWHSTL,RDRESL,LRRSOL,LRSYMI,LRSYMM C .. C .. Intrinsic Functions .. INTRINSIC ABS,NINT,REAL,SQRT C .. C .. Common blocks .. COMMON /ATSYM/RSYM(4,4,96),RSYMT(4,4,96),NSYM,PERM(4,4),IPRINT, + MATOMS,LATOMS,JX,JY,JZ,LPRINT,LASYM,NSG COMMON /CNTROL/NBIN,ST2MIN,DELST2,NPS,MIR2,IOPT,SMAX,IERR,WMAX, + WMAX2,SMIN,SMINCB,SMAXCB,IRFSIG COMMON /IO/ISYSR,ISYSW,NCOL,CONV,NMON COMMON /POSN/IFO,IFS,IAO,IW,ICOA,ICOB,ICOC,ICOD,IFC1,IAC1,IFC2, + IAC2,IFC3,IAC3,IAO2,IW2,ICOA2,ICOB2,ICOC2,ICOD2 COMMON /SIGS/SIGMAA(3,50),SIGMAN(3,50),SIGMAP(3,50),ACHANG(3) COMMON /TABS/NSTEP,COSA(120),SINA(120),ETAB(4001),ARGMAX,DLIM(51) C .. C .. Save statement .. SAVE C .. C .. Data statements .. C C NLPRGI gives number of LSPRGI C LSPRGI - Label String for PRoGram Input C CTPRGI - Column Type for PRoGram Input C C these columns may be added to original columns C NLPRGO gives number of LSPRGO C LSPRGO - Label String for PRoGram Output C CTPRGO - Column Type for PRoGram Output C DATA NLPRGI,LSPRGI/25,'H','K','L','FP','SIGFP','PHIBP','WP', + 'HLA','HLB','HLC','HLD','FC','PHIC','FC1','PHIC1', + 'FC2','PHIC2','FC3','PHIC3', + 'PHIB2','W2','HLA2','HLB2','HLC2','HLD2',475*' '/ DATA CTPRGI/'H','H','H','F','Q','P','W','A','A','A','A','F','P', + 'F','P','F','P','F','P','P','W','A','A','A','A',475*' '/ C DATA NLPRGO,LSPRGO/6,'WWT','DELFWT','FWT','PHWT','PHCMB','WCMB', C + 194*' '/ C DATA CTPRGO/'W','F','F','P','P','W',194*' '/ C If getting Sim weights for partial structure: output will be DELFWT...WCMB C If combining phases and partial structure: output will be DELFWT...WCMB PHCMB C If combining phases from 2 MIR sets of HL coeffs: output will be HLAC..HLDC WCMB PHCMB DATA NLPRGO,LSPRGO/10,'DELFWT','PHDELFWT','FWT','PHFWT', + 'WCMB','PHCMB', 'HLAC','HLBC','HLCC','HLDC', + 490*' '/ DATA CTPRGO/'F','P','F','P','W','P','A','A','A','A',490*' '/ C C---- Key words C DATA KEYWRD/'TITL','RESO','PART','COMB','SIGM','SYMM','RANG', + 'ERRO','LABI','LABO','END '/ C .. C LABOUT_SAVE = ' ' IOPT = 1 NBIN = 20 NPS = 1 IERR = 0 NMON = 1000 ISYSR = 5 ISYSW = 6 IFAIL = 0 SMIN = 0.0 SMAX = 0.5 RMAX = 0.0 SMINCB = -1.0 SMAXCB = -1.0 NSMP = 0 IFO=0 IFS=0 IAO=0 IW=0 ICOA=0 ICOB=0 ICOC=0 ICOD=0 IFC1=0 IAC1=0 IFC2=0 IAC2=0 IFC3=0 IAC3=0 IAO2=0 IW2=0 ICOA2=0 ICOB2=0 ICOC2=0 ICOD2=0 ICHK=0 IRFSIG=0 C These can be used to DAMP Sigmaa weight for Partial contributions ACHANG(1) = 1.0 ACHANG(2) = 1.0 ACHANG(3) = 1.0 ICHANG = 0 C C ************************************************ CALL CCPDPN(ISYSR,'DATA','READONLY','F',LDUM,IFAIL) CALL CCPDPN(ISYSW,'PRINTER','PRINTER','F',LDUM,IFAIL) C ************************************************ C SLIM = 0.25 10 CONTINUE LINE = ' ' NTOK = NPAR C C *********************************************************** CALL PARSER(KEY,LINE,IBEG,IEND,ITYP,FVALUE,CVALUE,IDEC,NTOK,LEND, + .TRUE.) C *********************************************************** C C---- End of file? C IF (LEND) GO TO 150 C C---- Loop over possible key-words C DO 20 I = 1,NKEYS IF (KEY.EQ.KEYWRD(I)) GO TO 30 20 CONTINUE C IF(KEY .EQ. 'ANAL') THEN CALL CCPERR(1,' *** PHISTATS should be used to analyse phases') ENDIF CALL CCPERR(1, '*** Error - unrecognized keyword: ' // + LINE(IBEG(1) :IEND(1))) C ******************** C 30 CONTINUE C C* DATA KEYWRD/'TITL','RESO','PART','COMB','SIGM','SYMM','RANG', C* + 'ERRO','LABI',LABO','END '/ C GO TO ( 40, 50, 60, 70, 80, 90, 110, + 120, 130, 140, 150) I C C---- TITLE C 40 TITLE = LINE(IBEG(2) :) GO TO 10 C C---- RESOlution C =========== C C Resolution limits - either 4sinsq/lamdasq or d C 50 CONTINUE ITOK = 2 C C ************************************************* CALL RDRESO(ITOK,ITYP,FVALUE,NTOK,RMIN,RMAX,SMIN,SMAX) C ************************************************* C WRITE (ISYSW,FMT=6002) RMIN,RMAX,SMIN,SMAX 6002 FORMAT (//' Resolution limits in As ',2F10.2,/11X,'as 4sinsq/', + 'lsq ',2F10.5,/) C GO TO 10 C C PARTial C ========= C C---- Set flag for no combination C C Number of partial structures = 1 C 60 CONTINUE NPS = 1 IOPT = 1 ICOMB = 0 ICHANG = 0 C DAMP D1 ( D2 ( D3) ) N = 2 CALL CCPUPC(CVALUE(N)) IF (CVALUE(N) .EQ. 'DAMP') THEN ICHANG = 1 N = N+1 C ****************************** CALL GTPREA(N,ACHANG(1),NTOK,ITYP,FVALUE) C ****************************** ACHANG(2) = ACHANG(1) ACHANG(3) = ACHANG(2) C Is there a 2nd DAMP factor? IF(N .GE. NTOK) GO TO 61 N = N+1 C ****************************** CALL GTPREA(N,ACHANG(2),NTOK,ITYP,FVALUE) C ****************************** ACHANG(3) = ACHANG(2) C And a 3rd? IF(N .GE. NTOK) GO TO 61 N = N+1 C ****************************** CALL GTPREA(N,ACHANG(3),NTOK,ITYP,FVALUE) C ****************************** END IF 61 CONTINUE WRITE (ISYSW,FMT=6004) 6004 FORMAT (/' Calculation of Fourier coefficients from partial stru', + 'cture',/) IF(ICHANG.NE.0) + WRITE(6,'(/'' Partial weights will be multiplied by'',3F8.3)') + ACHANG GO TO 10 C C COMBine [PART ] [MIR2] [RESO ] C [DAMP [ [ ]]] C ======= C C---- Set flag to combine coefficients, optionally set number of partial C structures and resolution limits for phase combination C Note that, if present, NPS must precede resolution limits C 70 IOPT = 3 ICOMB = 1 C C---- Read number of partial structures C Number of partial structures defaults to 1 C NPS = 1 MIR2 = 0 C N = 1 71 N = N+1 IF (N .LE. NTOK) THEN IF (ITYP(N) .NE. 1) CALL CCPERR(1,' Error on COMBINE LINE') C String read CALL CCPUPC(CVALUE(N)) C RESO IF (CVALUE(N) .EQ. 'RESO') THEN ITOK = N+1 C Read resolution limits. Note: must supply both limits C unless end of line. C MDW: have replaced RDRESL with RDRESO because (a) former is C a pain with a variable number of tokens, and (b) latter C is more standard for coders and users. CALL RDRESO(ITOK,ITYP,FVALUE,NTOK,RMINCB, + RMAXCB,SMINCB,SMAXCB) C Treat single limit after 'RESOLUTION' as low resolution IF (SMINCB .LT. 0.0) THEN SMINCB = SMAXCB RMINCB = RMAXCB SMAXCB = -1.0 RMAXCB = -1.0 ENDIF N = N+2 GO TO 71 ENDIF C---- Read number of partial structures C PART IF (CVALUE(N) .EQ. 'PART') THEN ICHANG = 0 N = N+1 C ****************************** CALL GTPINT(N,NPS,NTOK,ITYP,FVALUE) C ****************************** ICHK = 1 GO TO 71 ENDIF C DAMP [ [ ] ] IF (CVALUE(N) .EQ. 'DAMP') THEN ICHANG = 1 N = N+1 C ****************************** CALL GTPREA(N,ACHANG(1),NTOK,ITYP,FVALUE) C ****************************** ACHANG(2) = ACHANG(1) ACHANG(3) = ACHANG(2) C Is there a 2nd DAMP factor? IF (N+1.GT.NTOK .OR. ITYP(N+1).EQ.1) GO TO 71 N = N+1 C ****************************** CALL GTPREA(N,ACHANG(2),NTOK,ITYP,FVALUE) C ****************************** ACHANG(3) = ACHANG(2) C And a 3rd? IF (N+1.GT.NTOK .OR. ITYP(N+1).EQ.1) GO TO 71 N = N+1 C ****************************** CALL GTPREA(N,ACHANG(3),NTOK,ITYP,FVALUE) C ****************************** GO TO 71 ENDIF C MIR2 IF (CVALUE(N) .EQ. 'MIR2') THEN IF(ICHK.EQ.0) NPS = 0 MIR2 = 1 GO TO 71 ENDIF ENDIF C C---- NPS is number of partial structures to be combined. C IF (NPS.GT.1) WRITE (ISYSW,FMT=6006) NPS 6006 FORMAT (' Number of partial structures =',I3) C IF(ICHANG.EQ.1) + WRITE(6,'(/'' Partial weights will be multiplied by'',3F8.3)') + ACHANG C IF ((NPS.LT.0) .OR. (NPS.GT.3)) CALL CCPERR(1, + 'Number of partial structures must be between 0 and 3') C IF (MIR2.GT.0) WRITE (ISYSW,FMT=6007) 6007 FORMAT (' Combination of 2 experimental phase sets') C WRITE (ISYSW,FMT=6008) 6008 FORMAT (/' Phase combination and Fourier coefficients',/) GO TO 10 C C Sigmaa values from elsewhere. ( comfirm NPSA and NBINA here) C ======================================= C 80 CONTINUE NPSA=NPS NBINA = NBIN ITOK = 2 IF(ITOK.LE.NTOK)CALL GTPINT(ITOK,NPSA,NTOK,ITYP,FVALUE) ITOK = 3 IF(ITOK.LE.NTOK)CALL GTPINT(ITOK,NBINA,NTOK,ITYP,FVALUE) WRITE(ISYSW,6010) (IPSA,IPSA=1,NPSA) 6010 FORMAT(/,' INPUT SIGMAA VALUES',/, $ ' BIN_NO PARTIAL STRUCTURE',/, 7X,3(I5,5X)) DO 82 IBINA=1,NBINA READ(ISYSR,*) ( SIGMAA(IPSA,IBINA),IPSA=1,NPSA) WRITE(ISYSW,6011) $ IBINA,(SIGMAA(IPSA,IBINA),IPSA=1,NPS) 6011 FORMAT(I7,3F10.5) 82 CONTINUE IRFSIG = 1 GO TO 10 C C---- SYMMetry NspaceGroup OR NAME OF SPACEGROUP C ============================================= C 90 CONTINUE ITOK = 2 C C ************************************************************ CALL RDSYMM(ITOK,LINE,IBEG,IEND,ITYP,FVALUE,NTOK,NAMSPG,NSPGRP, + PGNAME,NSYM,NSMP,RSYM) C ************************************************************ C GO TO 10 C C RANGes read number of ranges (bins), monitor interval C ====================================================== C C---- Monitoring NUMBER _ program prints lots of information C about every nmon th reflection C 110 CONTINUE C C ******************************* CALL GTPINT(2,NBIN,NTOK,ITYP,FVALUE) C ******************************* C IF (NBIN.GT.50) CALL CCPERR(1, + 'The number of resolution ranges must be <= 50') C C ******************************* CALL GTPINT(3,NMON,NTOK,ITYP,FVALUE) C ******************************* C WRITE (ISYSW,FMT=6012) NMON 6012 FORMAT (/' Program prints lots of information about every nmon', + ' th reflection - nmon = ',I6,/) GO TO 10 C C ERROR C ===== C C---- Set error estimation flag C 120 IERR = 1 GO TO 10 C C LABI command, read input labels C =============================== C 130 CONTINUE C C **************************************** CALL LKYIN(1,LSPRGI,NLPRGI,NTOK,LINE,IBEG,IEND) C **************************************** C GO TO 10 C C---- LABO command, read output labels C ================================ C 140 CONTINUE LABOUT_SAVE = LINE GO TO 10 C C End of input, END command read C 150 CONTINUE C Move this here - NBIN may be read AFTER reso limits.. ST2MIN = SMIN/4.0 DELST2 = (SMAX/4.0-ST2MIN)/NBIN C C IF (IOPT.EQ.1) WRITE (ISYSW,FMT=6016) C 6016 FORMAT (' Option PARTIAL -- Map coefficients (a) difference map ', C + 'mFo-Fc [ DELFWT exp(PHWT) ]',/' ', C + ' (b) 2mFo-D*Fc [ FWT exp(PHWT) ]',/' HK', C + 'LOUT contains columns H K L S ... DELFWT FWT PHWT') C IF (IOPT.EQ.3) WRITE (ISYSW,FMT=6018) C 6018 FORMAT (' Option COMBINE -- ',/' phase combination and co', C + 'mbined map',/' coefficients 2mFo-D*Fc [ FWT exp(', C + 'PHWT) ]',/' HKLOUT contains columns H K L S ... P', C + 'HCMB WWT FWT PHWT',/' If no Hendrickson-Lattman coeffici', C + 'ents in input, ',/' set C=D=0 A=Fobs*W*cos(Phmir) B=', C + '..',/) call ccp4h_summary_beg() IF (IOPT.EQ.1) WRITE (ISYSW,FMT=6016) 6016 FORMAT (' Option PARTIAL -- ',/, + ' Map coefficients (a) difference map;Fo-DFc PHFo-DFc',/, + ' (b) 2mFo-D*Fc PH2mFo-D*Fc',/, + ' HKLOUT contains columns H K L S ... ', + ' MFODC 2MFODFC WCMB=(SimWt)') C IF (IOPT.EQ.3 .AND. NPS.NE.0) WRITE (ISYSW,FMT=6018) 6018 FORMAT (' Option COMBINE with partial structure-- ',/, + ' Phase combination of several sets of information:',/, + ' Combined map coefficients ',/, + ' Map coefficients (a) difference map;Fo-DFc PHFo-DFc',/, + ' (b) 2mFo-D*Fc PH2mFo-D*Fc',/, + ' HKLOUT contains columns H K L S ... ', + ' MFODC PHMFODC 2MFODFC PH2MFODFC WCMB PHCMB') C IF (IOPT.EQ.3 .AND. NPS .EQ. 0 .AND. MIR2.GT.0) + WRITE (ISYSW,FMT=6019) 6019 FORMAT (' Option Combine two MIR sets-- ',/, + ' HKLOUT contains columns H K L S ... ', + ' HLAC HLBC HLCC HLDC WCMB PHCMB') C C---- Unless option 3 is chosen, nps must be 1 C IF (IOPT.LT.3) NPS = 1 C IF (IOPT.EQ.3 .AND. NPS.GT.0) WRITE (ISYSW,FMT=6022) NPS 6022 FORMAT (' SigmaA determination for ',I2,' partial structure(s)') call ccp4h_summary_end() C C---- Set lookup to search out compulsory and optional columns C C Note LOOKUP = 0 optional, -1 compulsory C DO 170 I = 6,100 LOOKUP(I) = 0 170 CONTINUE C H K L FP SIGFP DO 180 I = 1,5 LOOKUP(I) = -1 180 CONTINUE C IF (IOPT.LT.3) THEN C C---- PARTial option, assign FC, PHIC C LOOKUP(12) = -1 LOOKUP(13) = -1 END IF C IF (IOPT.EQ.3) THEN C C COMBine option C C---- PHIB W for isomorphous phase, complusory C LOOKUP(6) = -1 LOOKUP(7) = -1 C C---- a b c d Hendrickson-Lattmann coefficients C C---- Not compulsory (soft) - set to c=d=0 a=eobs*W*cos(amir) b=') C LOOKUP(8) = 1 LOOKUP(9) = 1 LOOKUP(10) = 1 LOOKUP(11) = 1 C C---- FC PHIC is one partial structure C IF (NPS.GE.1) THEN C C FC PHIC C LOOKUP(12) = -1 LOOKUP(13) = -1 IF (NPS.GE.2) THEN C C FC1 PHCAL1 and FC2 PHCAL2 remove need for FC PHIC C LOOKUP(12) = 0 LOOKUP(13) = 0 LOOKUP(14) = -1 LOOKUP(15) = -1 LOOKUP(16) = -1 LOOKUP(17) = -1 IF (NPS.EQ.3) THEN C C FC3 PHCAL3 C LOOKUP(18) = -1 LOOKUP(19) = -1 END IF END IF END IF C C Do not set HLA2 ... as compulsory. IF (MIR2.EQ.1) THEN LOOKUP(20) = -1 LOOKUP(21) = -1 LOOKUP(22) = 1 LOOKUP(23) = 1 END IF END IF C C---- Print header info C MPRINT = 2 C C---- Open MTZ input file, return cell parameters, number of columns C and assign column labels. C C ************************************** CALL LROPEN(1,'HKLIN',MPRINT,IFAIL) CALL LRCELL(1,CELL) CALL LRCLAB(1,CLABS,CTYPS,NCOL) CALL LRASSN(1,LSPRGI,NLPRGI,LOOKUP,CTPRGI) C ************************************** IF (RMAX.EQ.0.0) THEN CALL LRRSOL(1,SMIN,SMAX) ST2MIN = SMIN/4.0 DELST2 = (SMAX/4.0-ST2MIN)/NBIN ENDIF C ST2MIN = SMIN/4.0 DELST2 = (SMAX/4.0-ST2MIN)/NBIN NBIN1 = NBIN + 1 C C Default resolution limits for phase combination equal to overall limits IF (SMINCB .LT. 0.0) THEN SMINCB = SMIN ENDIF IF (SMAXCB .LT. 0.0) THEN SMAXCB = SMAX ENDIF C DO 160 I = 1,NBIN1 STOL2 = REAL(I-1)*DELST2 + ST2MIN IF (STOL2.GT.0.0) DLIM(I) = 0.5/SQRT(STOL2) IF (STOL2.LE.0.0 .OR. DLIM(I).GT.1000.0) DLIM(I) = 999.99 160 CONTINUE C C---- To get back to Randy's form C STMIN = SQRT(ABS(ST2MIN)) STMAX = SQRT(ABS(SMAX/4.0)) WRITE (ISYSW,FMT=6024) STMIN,STMAX,NBIN 6024 FORMAT (/' Data from sin(theta)/lambda of',F8.4,' to',F8.4, + ' divided into',I4,' ranges') C IF (IOPT .EQ. 3) THEN IF (SMINCB .LT. SMIN) THEN WRITE (ISYSW, '(/A,A/)') $ ' **** WARNING: low resolution limit for', $ ' phase combination reset to overall limit ****' SMINCB = MAX(SMINCB, SMIN) ENDIF IF (SMAXCB .GT. SMAX) THEN WRITE (ISYSW, '(/A,A/)') $ ' **** WARNING: high resolution limit for', $ ' phase combination reset to overall limit ****' SMAXCB = MIN(SMAXCB, SMAX) ENDIF RMINCB = SQRT(1.0/SMINCB) RMAXCB = SQRT(1.0/SMAXCB) WRITE (ISYSW, '(/A,2F10.2,A/)') $ ' Resolution limits for phase combination are', $ RMINCB, RMAXCB, ' Angstrom' ENDIF C C ********************************************** CALL LRSYMI(1,NSYMP,LTYPE,NSPGRP,NAMSPG,PGNAME) IF(NSMP.GT.0 .AND.NSYMP.NE.NSMP) CALL CCPERR + (1,' Disagreement between symmetry operators') CALL LRSYMM(1,NSYM,RSYM) C ********************************************** NMULT = NSYM/NSYMP WRITE (ISYSW,FMT=6014) NMULT 6014 FORMAT (/' Lattice multiplicity',I6) C C---- Now I can work out centricity and epsilon classes C c DO 100 J = 1,NSYM C C ******************************** c CALL INVSYM(RSYM(1,1,J),RSYMT(1,1,J)) C ******************************** C c 100 CONTINUE C C ****************************** CALL CENTRIC(NSYM,RSYM,IPRINT) CALL EPSLN(NSYM,NSYMP,RSYM,IPRINT) C ****************************** C IFO = LOOKUP(4) IFS = LOOKUP(5) IAO = LOOKUP(6) IW = LOOKUP(7) ICOA = LOOKUP(8) ICOB = LOOKUP(9) ICOC = LOOKUP(10) ICOD = LOOKUP(11) IFC1 = LOOKUP(12) IAC1 = LOOKUP(13) C IF (NPS.GT.1) THEN IFC1 = LOOKUP(14) IAC1 = LOOKUP(15) IFC2 = LOOKUP(16) IAC2 = LOOKUP(17) IFC3 = LOOKUP(18) IAC3 = LOOKUP(19) END IF C C---- Add phase analysis and combination columns - C bypass Sim weighting if C these labels set. C IAO2 = LOOKUP(20) IW2 = LOOKUP(21) ICOA2 = LOOKUP(22) ICOB2 = LOOKUP(23) ICOC2 = LOOKUP(24) ICOD2 = LOOKUP(25) C C---- Leave cell dimensions unchanged C CELL(1) = 0.0 C C---- Read labels for output columns to add to h k l s C If iopt lt 3 output h k l FO ... plus Diff map "F" ,"2fo-dfc" , and W C phase same as PHIC, "F" and "2fo-dfc" may be + or - C C If iopt eq 3 and partial structure present C output h k l FO ... plus Diff map "F" and phase, C "2fo-dfc" and phase, and W, and comb phase C C If iopt eq 3 and no partial structure present ( ie two MIR sets) C output h k l FO ... plus HLAcomb HLBcomb .. W and comb phase C C NLPRGO = 4 C C IF (IOPT.EQ.3) THEN C LSPRGO(1) = LSPRGO(5) C LSPRGO(2) = LSPRGO(6) C CTPRGO(1) = CTPRGO(5) C CTPRGO(2) = CTPRGO(6) C END IF C C IF(IOPT.LT.3) NLPRGO = 5 IF(IOPT.LT.3) THEN C Partial structure - only need to output mFODFC and 2mFOdFc and weight; C - phase is always PHIC NLPRGO = 3 LSPRGO(2) = LSPRGO(3) LSPRGO(3) = LSPRGO(5) CTPRGO(2) = CTPRGO(3) CTPRGO(3) = CTPRGO(5) END IF IF(IOPT.EQ.3 .AND. NPS.NE.0 .AND. MIR2.EQ.0 ) NLPRGO = 6 IF(IOPT.EQ.3 .AND. NPS.NE.0 .AND. MIR2.GT.0 ) NLPRGO = 10 C Dont need partial map coeffs with bias removed if there is no partial structure IF(IOPT.EQ.3 .AND. NPS.EQ.0 .AND. MIR2.GT.0 ) THEN NLPRGO = 6 LSPRGO(1) = LSPRGO(7) LSPRGO(2) = LSPRGO(8) LSPRGO(3) = LSPRGO(9) LSPRGO(4) = LSPRGO(10) CTPRGO(1) = CTPRGO(7) CTPRGO(2) = CTPRGO(8) CTPRGO(3) = CTPRGO(9) CTPRGO(4) = CTPRGO(10) END IF C C If we are combining 2 MR sets output HLAcomb etc... CALL LWHSTL (1, ' ') C LINE = LABOUT_SAVE C C ****************************************************** CALL PARSER(KEY,LINE,IBEG,IEND,ITYP,FVALUE,CVALUE,IDEC,NTOK, + LEND,.TRUE.) CALL LKYOUT(1,LSPRGO,NLPRGO,NTOK,LINE,IBEG,IEND) C ****************************************************** C C---- Add to title and labels C C ******************************** CALL LWOPEN(1,'HKLOUT') CALL LWTITL(1,TITLE,1) CALL LWASSN(1,LSPRGO,NLPRGO,CTPRGO,1) C ******************************** END C ===================================== SUBROUTINE CENPHA(A,B,ALF,W) C ===================================== C ALF derived from A and B; ALF units are radians. C .. Scalar Arguments .. REAL A,ALF,B,W C .. C .. Local Scalars .. REAL TWOPI C .. C .. Intrinsic Functions .. INTRINSIC ATAN2,SQRT,TANH C .. C .. Save statement .. SAVE C .. DATA TWOPI/6.283185/ C .. Data statements .. C .. C C---- Return centric phase and figure of merit from hendrickson- C lattman a and b. C the argument for tanh corresponds to x/2 (cf. s/r srin) C W = TANH(SQRT(A**2+B**2)) C IF (W.EQ.0.0) ALF = 0.0 IF (W.NE.0.0) THEN ALF = ATAN2(B,A) IF (ALF.LT.0.0) ALF = ALF + TWOPI C Cdebug WRITE(6,9876)A,B,C,D,ALP,W C9876 FORMAT(' IN CENTPH A B C D ALF W',4F8.3,5X,2F8.3) C END IF END C ================= SUBROUTINE COMBIN C ================= C C combine phases (0 to 3 partial structures and 0 or 1 C set of hendrickson-lattman coefficients for non-structural C phase information), compute map coefficients to minimize C model bias (explained in s/r unbias), and generate an even C more excessive quantity of statistics. C C---- Variables for combined phase stats. subdivided to partial C structure, mir and combined phase for centric and non-centric C data. only includes those reflections for which both C sources are combined. variation of information, cosine of C phase difference and partial structure weight stuff only C tabulated for non-centric data. C C---- Variables for stats on reflections using only partial C structure phases. centric and non-centric data. C C---- Variables for stats on reflections using only mir phases. C centric and non-centric data. C C .. Parameters .. INTEGER MCOLS PARAMETER (MCOLS=500) C .. C .. Scalars in Common .. REAL ARGMAX,CONV,DELST2,SMAX,SMAXCB,SMIN,SMINCB,ST2MIN,WMAX,WMAX2 INTEGER IAC1,IAC2,IAC3,IAO,IAO2,ICOA,ICOA2,ICOB,ICOB2,ICOC,ICOC2, + ICOD,ICOD2,IERR,IFC1,IFC2,IFC3,IFO,IFS,IOPT,ISYSR,ISYSW, + IW,IW2,NBIN,NCOL,NMON,NPS,NSTEP,MIR2,IRFSIG C .. C .. Arrays in Common .. REAL COSA,DLIM,ETAB, + SIGMAA,SIGMAN,SIGMAP,SINA,ACHANG C .. C .. Local Scalars .. REAL A,ALFCOM,ALFOUT,ALFOUTD,AO,AO2,AOUT,B,BOUT,C,D,EO,EPSI,FO, + FOUT,OVERM,OVERMC,OVERMN,PHI,SSQ,STOL2,W,WCHK,WCOM,WIR,WIR2,X REAL A2,B2,C2,D2,AOUTD,BOUTD,FOUTD INTEGER I,I1,IAC,IC,ICHK,IFC,INDMIR,INDMR2,INDPAR,INDX,IPS,ISYSAB, + J,NC,NN,NREF,NSHL,NT LOGICAL EOF C .. C .. Local Arrays .. REAL ABTAB(21),ADATA(MCOLS),APAR(3),BPAR(3),EC(3),FC(3), + R(11) INTEGER HIST(10,50),HISTT(10),IH(3) LOGICAL LOGMSS(MCOLS) C .. C .. External Subroutines .. EXTERNAL CENTR,COMFAZ,EPSLON,LRREFL,LWREFL,MIRFAZ,PARFAZ, + STATC,STATM,STATP,STINIT C .. C .. Intrinsic Functions .. INTRINSIC ATAN2,COS,INT,LOG,MIN,MOD,NINT,REAL,SIN,SQRT C .. C .. Common blocks .. COMMON /CNTROL/NBIN,ST2MIN,DELST2,NPS,MIR2,IOPT,SMAX,IERR,WMAX, + WMAX2,SMIN,SMINCB,SMAXCB,IRFSIG COMMON /IO/ISYSR,ISYSW,NCOL,CONV,NMON COMMON /POSN/IFO,IFS,IAO,IW,ICOA,ICOB,ICOC,ICOD,IFC1,IAC1,IFC2, + IAC2,IFC3,IAC3,IAO2,IW2,ICOA2,ICOB2,ICOC2,ICOD2 COMMON /SIGS/SIGMAA(3,50),SIGMAN(3,50),SIGMAP(3,50),ACHANG(3) COMMON /TABS/NSTEP,COSA(120),SINA(120),ETAB(4001),ARGMAX,DLIM(51) C .. C .. Save statement .. SAVE C .. C C Nicked from COMBINE C Set up table for getting A and B from W and phase. C Set A = X*cos(amir) B = X*sin(amir) C where W = I1(X)/I0(X) C X IN RANGE 0's) C c CW = RDUMP(2) c INDX = INT(CW*10.0) + 1 INDX = INT(WCOM*10.0) + 1 IF (INDX.GT.10) INDX = 10 HIST(INDX,NSHL) = HIST(INDX,NSHL) + 1 HISTT(INDX) = HISTT(INDX) + 1 c IF (IC.EQ.0) OVERMN = OVERMN + CW c IF (IC.EQ.0) NN = NN + 1 c IF (IC.NE.0) OVERMC = OVERMC + CW IF (IC.EQ.0) OVERMN = OVERMN + WCOM IF (IC.EQ.0) NN = NN + 1 IF (IC.NE.0) OVERMC = OVERMC + WCOM IF (IC.NE.0) NC = NC + 1 C WCOM = 0 AOUT = 0 BOUT = 0 FOUT = 0 AOUTD = 0 BOUTD = 0 FOUTD = 0 INDPAR = 0 INDMIR = 0 INDMR2 = 0 GO TO 20 END IF C C---- Put out some diagnostics. first partial structure phase only, C mir phase only, then combined phase stats. C 80 CONTINUE call ccp4h_summary_beg() WRITE (ISYSW,FMT=6006) 6006 FORMAT (' Calculation of Map Coefficients is COMPLETE') call ccp4h_summary_end() WRITE (ISYSW,FMT=6007) 6007 FORMAT (' Statistics on Phases and Figures of MERIT:') C C ******************** CALL STATP(NPS,NBIN,DLIM) CALL STATM(NBIN,DLIM) CALL STATC(NPS,NBIN,DLIM) C ******************** C C---- Now overall stats C DO 90 I = 1,11 R(I) = REAL(I-1)/10.0 90 CONTINUE C WRITE (ISYSW,FMT=6008) R 6008 FORMAT (' ',/' Figure of Merit Histograms (all DATA)',/' D L', + 'imits',F4.1,10F5.1) C DO 100 I = 1,NBIN I1 = I + 1 WRITE (ISYSW,FMT=6010) DLIM(I),DLIM(I1), (HIST(J,I),J=1,10) 6010 FORMAT (F7.2,'--',F5.2,10I5) 100 CONTINUE C WRITE (ISYSW,FMT=6012) HISTT 6012 FORMAT (' Overall ',10I5) NT = NN + NC IF (NT.NE.0) OVERM = (OVERMN+OVERMC)/NT IF (NN.NE.0) OVERMN = OVERMN/NN IF (NC.NE.0) OVERMC = OVERMC/NC WRITE (ISYSW,FMT=6014) NC,OVERMC,NN,OVERMN,NT,OVERM 6014 FORMAT (' ',/' Overall Figures of Merit',/' DATA N ', + '',/' CENTRIC ',I6,F10.5,/' NON-CENTRIC',I6,F10.5, + /' TOTAL ',I6,F10.5) RETURN C END C =============================================================== SUBROUTINE COMFAZ(IH,NPS,IC,NSHL,FO,FC,APAR,BPAR,A,B,C,D,ALFCOM, + WCOM,AOUT,BOUT,AOUTD,BOUTD) C =============================================================== C C C .. Scalar Arguments .. REAL A,AOUT,AOUTD,B,BOUT,BOUTD,C,D,FO,ALFCOM,WCOM INTEGER IC,NPS,NSHL C .. C .. Array Arguments .. REAL ALFPAR(3),APAR(3),BPAR(3),FC(3),WPAR(3),VARP(3) INTEGER IH(3) C .. C .. Arrays in Common .. REAL CDMIR,CDPC,SMMCOM,SMCMC,SMCPC,SMNCOM,SMNMC,SMNPC, + SIGMAA,SIGMAN,SIGMAP,VARPC,VARCOM,VARMIR,WEIGHT,ACHANG INTEGER NCCOM,NNCOM C .. C .. Local Scalars .. REAL ALF,ALFMIR,ATOT,BTOT,DLUZ,VARINF,VARM,W,WTDEN C .. C .. External Subroutines .. EXTERNAL CENPHA,MIRFIL,PHASE C .. C .. Intrinsic Functions .. INTRINSIC COS,SIN C .. C .. Common blocks .. COMMON /SIGS/SIGMAA(3,50),SIGMAN(3,50),SIGMAP(3,50),ACHANG(3) COMMON /STCOM/NCCOM(50),NNCOM(50),SMCPC(3,50),SMNPC(3,50), + VARPC(3,50),CDPC(3,50),SMCMC(50),SMNMC(50),VARMIR(50), + CDMIR(50),SMMCOM(50),SMNCOM(50),VARCOM(50),WEIGHT(3,50) C .. C .. Save statement .. SAVE C .. C IF (IC.EQ.0) NNCOM(NSHL) = NNCOM(NSHL) + 1 IF (IC.NE.0) NCCOM(NSHL) = NCCOM(NSHL) + 1 ATOT = A BTOT = B WTDEN = 0.0 C C---- Fill in rdump with partial structure and mir phase info C C ******************************************************** c CALL PARFIL(IH,NPS,IC,NSHL,FC,APAR,BPAR,RDUMP,ATOT,BTOT,WTDEN, c + SMCPC,SMNPC,VARPC) CALL PARFIL(IH,NPS,IC,NSHL,APAR,BPAR,ALFPAR,WPAR,VARP, + ATOT,BTOT,WTDEN,SMCPC,SMNPC,VARPC) c CALL MIRFIL(IH,NPS,IC,NSHL,A,B,C,D,RDUMP,SMCMC,SMNMC) CALL MIRFIL(IH,NPS,IC,NSHL,A,B,C,D,ALF,W,VARM,SMCMC,SMNMC) C ******************************************************** C c MIROFF = 4*NPS + 3 c VARM = RDUMP(MIROFF+7) WTDEN = WTDEN + VARM ALFMIR = ALF C Cdebug write(6,9876)ih,ATOT,BTOT,C,D C9876 format(' call phase-comfaz- h k l a b c d',3i3,4f8.3,5x,2f8.3) C C ***************************************** IF (IC.EQ.0) CALL PHASE(IH,ATOT,BTOT,C,D,ALFCOM,WCOM,VARINF) IF (IC.NE.0) CALL CENPHA(ATOT,BTOT,ALFCOM,WCOM) C ***************************************** C c RDUMP(1) = ALF c RDUMP(2) = W C IF (IC.EQ.0) RDUMP(3) = VARINF AOUT = WCOM*FO*COS(ALFCOM) BOUT = WCOM*FO*SIN(ALFCOM) C C---- Remove partial structure bias from map. - Randys comments C NOW WE HAVE COMBINED PHASE WEIGHTED FO. THIS IS ALL WE NEED C FOR CENTRIC OPTION 3 COEFFICIENTS. FOR NON-CENTRIC OPTION 3, C REMOVE PARTIAL STRUCTURE BIAS WITH UNBIAS. FOR NON-CENTRIC C OPTION 4, UNBIAS JUST KEEPS SOME STATISTICS. FOR OPTION 4, C TAKE THE DIFFERENCE WITH PARTIAL STRUCTURE 1 HERE. C DLUZ = SIGMAA(1,NSHL) * SQRT(SIGMAN(1,NSHL)/SIGMAP(1,NSHL)) AOUTD = AOUT - DLUZ*FC(1)*COS(ALFPAR(1)) BOUTD = BOUT - DLUZ*FC(1)*SIN(ALFPAR(1)) C C ************************************************* c IF (IC.EQ.0) CALL UNBIAS(ALF,RDUMP,WTDEN,NPS,NSHL,AOUT,BOUT, c + WEIGHT,CDPC) c IF (IC.EQ.0) CALL UNBIAS(ALFCOM,FC,ALFP,WP,WTDEN,NPS,NSHL,AOUT, CALL UNBIAS(IC,ALFCOM,FC,ALFPAR,VARP,WTDEN,NPS,NSHL, + AOUT,BOUT,WEIGHT,CDPC) C C ************************************************* C C---- Accumulate some sums, mostly for non-centric data C IF (IC.NE.0) SMMCOM(NSHL) = SMMCOM(NSHL) + WCOM C IF (IC.EQ.0) THEN SMNCOM(NSHL) = SMNCOM(NSHL) + WCOM VARCOM(NSHL) = VARCOM(NSHL) + VARINF VARMIR(NSHL) = VARMIR(NSHL) + VARM C ALFM = RDUMP(MIROFF+5) CDMIR(NSHL) = COS(ALFMIR-ALFCOM) + CDMIR(NSHL) END IF END C =============================================================== SUBROUTINE COMFZ2(IH,NPS,IC,NSHL,FO,A2,B2,C2,D2,A,B,C,D,ALFCOM, + WCOM) C =============================================================== C C C .. Scalar Arguments .. REAL A,A2,B,B2,C,C2,D,D2,FO,ALFCOM,WCOM INTEGER IC,NPS,NSHL C .. C .. Array Arguments .. INTEGER IH(3) C .. C .. Arrays in Common .. REAL CDMIR,CDPC,SMMCOM,SMCMC,SMCPC,SMNCOM,SMNMC,SMNPC,VARCOM, + VARMIR,VARPC,WEIGHT INTEGER NCCOM,NNCOM C .. C .. Local Scalars .. REAL ALF1,ATOT,BTOT,CTOT,DTOT,VARINF,VARM1,W1,WTDEN C .. C .. External Subroutines .. EXTERNAL CENPHA,MIRFIL,PARFIL,PHASE,UNBIAS C .. C .. Intrinsic Functions .. INTRINSIC COS,SIN C .. C .. Common blocks .. COMMON /STCOM/NCCOM(50),NNCOM(50),SMCPC(3,50),SMNPC(3,50), + VARPC(3,50),CDPC(3,50),SMCMC(50),SMNMC(50),VARMIR(50), + CDMIR(50),SMMCOM(50),SMNCOM(50),VARCOM(50),WEIGHT(3,50) C .. C .. Save statement .. SAVE C .. C IF (IC.EQ.0) NNCOM(NSHL) = NNCOM(NSHL) + 1 IF (IC.NE.0) NCCOM(NSHL) = NCCOM(NSHL) + 1 CALL MIRFIL(IH,NPS,IC,NSHL,A,B,C,D,ALF1,W1,VARM1,SMCMC,SMNMC) ATOT = A + A2 BTOT = B + B2 CTOT = C + C2 DTOT = D + D2 WTDEN = 0.0 C C ******************************************************** C Cdebug write(6,9876)ih,ATOT,BTOT,C,D C9876 format(' call phase-comfaz- h k l a b c d',3i3,4f8.3,5x,2f8.3) C C ***************************************** IF (IC.EQ.0) CALL PHASE(IH,ATOT,BTOT,CTOT,DTOT,ALFCOM,WCOM,VARINF) IF (IC.NE.0) CALL CENPHA(ATOT,BTOT,ALFCOM,WCOM) C ***************************************** WTDEN = WTDEN + VARINF C c RDUMP(1) = ALF c RDUMP(2) = W C IF (IC.EQ.0) RDUMP(3) = VARINF C C---- Accumulate some sums, mostly for non-centric data C IF (IC.NE.0) SMMCOM(NSHL) = SMMCOM(NSHL) + WCOM C IF (IC.EQ.0) THEN SMNCOM(NSHL) = SMNCOM(NSHL) + WCOM VARCOM(NSHL) = VARCOM(NSHL) + VARINF VARMIR(NSHL) = VARMIR(NSHL) + VARM1 C ALFM = RDUMP(MIROFF+5) CDMIR(NSHL) = COS(ALF1-ALFCOM) + CDMIR(NSHL) END IF END C ======================= SUBROUTINE DETSIG(IPS0) C ======================= C C .. Parameters .. INTEGER ISIZ PARAMETER (ISIZ=1000000) INTEGER MCOLS PARAMETER (MCOLS=500) C .. C .. Scalar Arguments .. INTEGER IPS0 C .. C .. Scalars in Common .. REAL ARGMAX,CONV,DELST2,SMAX,SMAXCB,SMIN,SMINCB,ST2MIN,WMAX,WMAX2 INTEGER IAC1,IAC2,IAC3,IAO,IAO2,ICOA,ICOA2,ICOB,ICOB2,ICOC,ICOC2, + ICOD,ICOD2,IERR,IFC1,IFC2,IFC3,IFO,IFS,IOPT,ISYSR,ISYSW, + IW,IW2,NBIN,NCOL,NMON,NPS,NSTEP,MIR2,IRFSIG C .. C .. Arrays in Common .. REAL COSA,DLIM,ETAB,SIGMAA,SIGMAN,SIGMAP,SINA,ACHANG C .. C .. Local Scalars .. REAL AVE02,AVE04,AVE20,AVE22,AVE40,CCDEN,CCNUM,CCTOT,DEL,DELMAX, + EPSI,FC,FO,OLD,OVERM,SIGA,SIGDEN,SIGLOG,SIGNUM,SSQ,STHOL, + STOL2,SUM02T,SUM04T,SUM20T,SUM22T,SUM40T,W,WDX,WT,X INTEGER I,I1,IC,ICHK,IFC,II,IPS,ISYSAB,NCHK,NRECOM, + NREF,NSHL,NTOT,NXY LOGICAL EOF,LNAN CHARACTER CXY*1 C .. C .. Local Arrays .. DOUBLE PRECISION DR(50),R(50) REAL ADATA(MCOLS),ASTL2(50),EC(ISIZ),EO(ISIZ),SUM02(50),SUM04(50), + SUM20(50),SUM22(50),SUM40(50),SUMM(50),SUMW(50),XY(2,50) INTEGER IBIN(ISIZ),IH(3),NE(50) CHARACTER LEGEND(20)*1,TITLEX(20)*1,TITLEY(20)*1 LOGICAL LOGMSS(MCOLS) C .. C .. External Subroutines .. EXTERNAL CCPERR,CENTR,EPSLON,LRREFL,PLOTR,SRIN C .. C .. Intrinsic Functions .. INTRINSIC ABS,DBLE,INT,LOG,MOD,NINT,REAL,SQRT C .. C .. Common blocks .. COMMON /CNTROL/NBIN,ST2MIN,DELST2,NPS,MIR2,IOPT,SMAX,IERR,WMAX, + WMAX2,SMIN,SMINCB,SMAXCB,IRFSIG COMMON /IO/ISYSR,ISYSW,NCOL,CONV,NMON COMMON /POSN/IFO,IFS,IAO,IW,ICOA,ICOB,ICOC,ICOD,IFC1,IAC1,IFC2, + IAC2,IFC3,IAC3,IAO2,IW2,ICOA2,ICOB2,ICOC2,ICOD2 COMMON /SIGS/SIGMAA(3,50),SIGMAN(3,50),SIGMAP(3,50),ACHANG(3) COMMON /TABS/NSTEP,COSA(120),SINA(120),ETAB(4001),ARGMAX,DLIM(51) C .. C .. Save statement .. SAVE C .. C .. Data statements .. DATA CXY/'*'/ DATA TITLEX/3*' ','(','S','I','N','(','T','H',')','/','L',')','*', + '*','2',3*' '/ DATA TITLEY/5*' ','L','N','(','S','I','G','M','A','A',')',5*' '/ DATA LEGEND/5*' ','S','I','G','M','A','A',' ','P','L','O','T', + 4*' '/ C .. LNAN = .TRUE. C C---- Find maxima for W and W2 C WMAX = -999999.0 WMAX2 = -999999.0 C C---- Deal with No partial structures case. C IPS = IPS0 IF (IPS.EQ.0) IPS = 1 SMAX = 4.0*DELST2*NBIN + 4.0*ST2MIN DO 10 I = 1,NBIN NE(I) = 0 SIGMAN(IPS,I) = 0.0000001 SIGMAP(IPS,I) = 0.0000001 SUM22(I) = 0.0 SUM20(I) = 0.0 SUM02(I) = 0.0 SUM40(I) = 0.0 SUM04(I) = 0.0 SUMW(I) = 0.0 ASTL2(I) = 0.0 10 CONTINUE NREF = 0 FC = 0.0 C C---- Read in data; store fo and fc (corrected for multiplicity but C not yet normalized) in eo and ec; store resolution bin number C in ibin (set negative as flag for centrics); accumulate sums C for normalization. weight |f|**2 for normalization by factor C of 1 for centrics, 2 for non-centrics. this makes the C similarly weighted sigmaa estimates come out more simply in C the refinement below, because the weighted mean |e|**2 is 1. C IF (IPS.EQ.1) IFC = IFC1 IF (IPS.EQ.2) IFC = IFC2 IF (IPS.EQ.3) IFC = IFC3 WRITE (ISYSW,FMT=6000) 6000 FORMAT (///) WRITE (ISYSW,FMT=*) + ' DETSIG-REF NO - h k l FO FC NSHL IC EPSI' 20 CONTINUE C C---- This reads one piece of data from the mtz file C and returns it in adata - (100 real array) C fp will return in adata(4),fc1 in adata(12),etc C C *********************** CALL LRREFL(1,SSQ,ADATA,EOF) C *********************** C IF (EOF) GO TO 30 CALL LRREFM(1,LOGMSS) C C---- Check sd(f) and resolution limits C IF (LOGMSS(IFS) .OR. LOGMSS(IFO)) THEN IF (LNAN) THEN IF (LOGMSS(IFO) .NEQV. LOGMSS(IFS)) + WRITE(ISYSW,FMT='(/,A,/,A,/)') + ' *** Warning: either FO/SIGFO=MNF but other is not.', + ' This must mean the column pair are inconsistent ***' LNAN = .FALSE. ENDIF GO TO 20 ELSE IF (ADATA(IFS).LT.0.000001) THEN GO TO 20 ELSE IF (SSQ.GT.SMAX .OR. SSQ.LT.SMIN) THEN GO TO 20 ELSE NREF = NREF + 1 IF (NREF.GT.ISIZ) THEN GO TO 160 ELSE C IH(1) = NINT(ADATA(1)) IH(2) = NINT(ADATA(2)) IH(3) = NINT(ADATA(3)) W = 0.0 IF (IW.NE.0) THEN IF (.NOT.LOGMSS(IW)) THEN IF (ADATA(IW).GT.0.0) W = ADATA(IW) ENDIF END IF IF (W.GT.WMAX) WMAX = W W = 0.0 IF (IW2.NE.0) THEN IF (.NOT.LOGMSS(IW2)) THEN IF (ADATA(IW2).GT.0.0) W = ADATA(IW2) ENDIF END IF IF (W.GT.WMAX2) WMAX2 = W C FO = ADATA(IFO) C Should bail out if no FC? IF (IFC.GT.0)THEN IF(LOGMSS(IFC)) THEN WRITE(ISYSW,FMT='(/,A,3I5,/)') + ' *** WARNING: this reflection has no FC ***',IH NREF = NREF - 1 GO TO 20 ENDIF FC = ADATA(IFC) END IF STOL2 = SSQ/4.0 NSHL = INT((STOL2-ST2MIN)/DELST2) + 1 IF (NSHL.LT.1) NSHL = 1 IF (NSHL.GT.NBIN) NSHL = NBIN C C ************ CALL CENTR(IH,IC) C ************ C IBIN(NREF) = NSHL IF (IC.EQ.1) IBIN(NREF) = -NSHL NE(NSHL) = NE(NSHL) + 1 ASTL2(NSHL) = ASTL2(NSHL) + STOL2 C C ********************** CALL EPSLON(IH,EPSI,ISYSAB) C ********************** C ICHK = MOD(NREF,NMON) IF (ICHK.EQ.1) WRITE (ISYSW,FMT=6002) IPS,NREF,IH,FO,FC,NSHL, + IC,EPSI 6002 FORMAT (2I7,3X,3I3,4X,2F6.1,2I4,F4.0) C WT = 2.0 - IC SUMW(NSHL) = SUMW(NSHL) + WT SIGMAN(IPS,NSHL) = FO**2*WT/EPSI + SIGMAN(IPS,NSHL) SIGMAP(IPS,NSHL) = FC**2*WT/EPSI + SIGMAP(IPS,NSHL) EO(NREF) = FO/SQRT(EPSI) EC(NREF) = FC/SQRT(EPSI) GO TO 20 END IF END IF C C---- End of input C 30 CONTINUE C C---- Take averages. WRITE (6,FMT='(/,a,f12.3,a,f12.3)') ' Maximum WP ',WMAX, + ' MAXIMUM W2 ',WMAX2 C NCHK = NREF/NBIN NRECOM = 500 IF (NCHK.LT.NRECOM) WRITE (6,FMT='(/,A,/,A,I8,/,A,I8,/)') + ' ****WARNING!!! - YOU HAVE TOO FEW REFLECTIONS PER BIN ****', + ' Average number per bin is ',NCHK, + ' Recommended number is at least ',NRECOM DO 40 I = 1,NBIN IF (NE(I).NE.0) THEN SIGMAN(IPS,I) = SIGMAN(IPS,I)/SUMW(I) SIGMAP(IPS,I) = SIGMAP(IPS,I)/SUMW(I) ASTL2(I) = ASTL2(I)/NE(I) END IF 40 CONTINUE IF (IPS0.NE.0) THEN C C---- Normalize to get e's and take sums for correlation coefficient C between eo**2 and ec**2. initial estimate for sigmaa is sqrt C of this cc (hauptman's alpha in acta cryst a38: 289-294(1982)) C also get overall correlation between e**2's. C DO 50 I = 1,NREF NSHL = ABS(IBIN(I)) ICHK = MOD(NREF,NMON) EO(I) = EO(I)/SQRT(SIGMAN(IPS,NSHL)) EC(I) = EC(I)/SQRT(SIGMAP(IPS,NSHL)) SUM22(NSHL) = (EO(I)*EC(I))**2 + SUM22(NSHL) SUM20(NSHL) = EO(I)**2 + SUM20(NSHL) SUM02(NSHL) = EC(I)**2 + SUM02(NSHL) SUM40(NSHL) = EO(I)**4 + SUM40(NSHL) SUM04(NSHL) = EC(I)**4 + SUM04(NSHL) 50 CONTINUE SUM22T = 0.0 SUM20T = 0.0 SUM02T = 0.0 SUM40T = 0.0 SUM04T = 0.0 NTOT = 0 IF ( IRFSIG.NE.0) WRITE(ISYSW,FMT=6003) 6003 FORMAT(/' ****************************************************',/, $ /' Input SIGMAA values will be used, instead of the',/, $ ' following initial estimates from the data statistics', $ /' ****************************************************') WRITE (ISYSW,FMT=6004) 6004 FORMAT (' Initial Estimates for SIGMAA',/' N STHOL D Li', + 'mits Sigman Sigmap <(EO*EC)**2> SIGMAA SIG_init') C DO 60 I = 1,NBIN C Use XY array to hold SIGMAA incase it is preset.. XY(1,I) = 0.0 IF ( IRFSIG.NE.0) XY(1,I) = SIGMAA(IPS,I) IF (NE(I).NE.0) THEN NTOT = NE(I) + NTOT SUM22T = SUM22(I) + SUM22T SUM20T = SUM20(I) + SUM20T SUM02T = SUM02(I) + SUM02T SUM40T = SUM40(I) + SUM40T SUM04T = SUM04(I) + SUM04T SIGNUM = NE(I)*SUM22(I) - SUM20(I)*SUM02(I) IF (SIGNUM.GT.0.0) THEN SIGDEN = SQRT((NE(I)*SUM40(I)-SUM20(I)**2)* + (NE(I)*SUM04(I)-SUM02(I)**2)) C C---- Signum/sigden is correlation between e**2's. the initial C estimate of sigmaa is the sqrt of this correlation C coefficient(hauptman, 1982). C SIGMAA(IPS,I) = SQRT(SIGNUM/SIGDEN) ELSE WRITE (ISYSW,FMT=6006) 6006 FORMAT (' Correlation between E**2''S is Non-Positive --', + ' SIGMAA set to 0.05') SIGMAA(IPS,I) = 0.05 END IF AVE20 = SUM20(I)/NE(I) AVE02 = SUM02(I)/NE(I) AVE40 = SUM40(I)/NE(I) AVE04 = SUM04(I)/NE(I) AVE22 = SUM22(I)/NE(I) STHOL = SQRT(ASTL2(I)) I1 = I + 1 WRITE (ISYSW,FMT=6008) NE(I),STHOL,DLIM(I),DLIM(I1), + SIGMAN(IPS,I),SIGMAP(IPS,I),AVE20,AVE02,AVE40,AVE04,AVE22, + SIGMAA(IPS,I),XY(1,I) 6008 FORMAT (I6,F7.4,F6.2,'--',F5.2,1P,2E11.4,0P,4F9.4,F12.5, + 2F9.4) END IF IF ( IRFSIG.NE.0) SIGMAA(IPS,I) = XY(1,I) 60 CONTINUE CCNUM = NTOT*SUM22T - SUM20T*SUM02T CCDEN = SQRT((NTOT*SUM40T-SUM20T**2)* (NTOT*SUM04T-SUM02T**2)) C PJX: There is a problem if for some reason CCDEN turns out to C be zero C This fix suggested by Dusan Turk - like him I haven't checked C whether setting CCTOT to zero in these circumstances is a C crystallographically sensible thing to do, although he does say C "after the check for the zero divide trap does give reasonable C results" IF (ABS(CCDEN).GT.0.001) THEN CCTOT = CCNUM/CCDEN ELSE CALL CCPERR(2, + 'S/R DETSIG: CCDEN < 0.001: setting overall CC to 0.0') WRITE (ISYSW,FMT='(7(/,A))') + ' In CCTOT calculation: CCNUM = ',CCNUM, + ' CCDEN = ',CCDEN, + ' NTOT = ',NTOT, + ' SUM22T= ',SUM22T, + ' SUM20T= ',SUM20T, + ' SUM40T= ',SUM40T, + ' SUM04T= ',SUM04T WRITE (ISYSW,FMT='(3(/,A))') + ' ** The Overall Correlation has been set to zero **', + ' ** You must check that the results of the sigmaA **', + ' ** procedure are still reasonable in this case! **' CCTOT = 0.0 END IF WRITE (ISYSW,FMT=6010) NTOT,CCTOT 6010 FORMAT (' Overall Correlation on |E|**2 for',I6,' Reflections ', + 'is',F10.5) C C If we are not determining sigmaa from the data, then we're done! C IF (IRFSIG.NE.0) RETURN C C---- Iteratively improve estimate of sigmaa. C sigmaa = <|eo|*|ec|*cos(da)> = <|eo|*|ec|*m> (srinivasan and C chandrasekaran, indian j. pure appl. phys. 4:178 (1966)), and C m is a function of sigmaa. with the correct sigmaa, the mean C value of (sigmaa-<|eo|*|ec|*m>) should be zero. with data C normalized so that weighted |e|**2 is equal to 1 and C weighting the centric(1) and non-centric(2) terms, this C corresponds to the maximum likelihood estimate of lunin and C urzhumtsev (acta cryst. a40:269 (1984)), although they C considered only non-centric reflections. use newton's C method to find zero. form of derivative is slightly different C for centric and non-centric because of different expressions C for m (see s/r srin). iterate until no shift is C greater than 0.0001, or for a maximum of 10 cycles. C DO 100 II = 1,10 WRITE (ISYSW,FMT=6012) II 6012 FORMAT (' Iteration Number',I3,' on SIGMAA',/' ', + ' OLD NEW MEAN',/' D Limi', + 'ts SIGMAA Shift SIGMAA W') C DO 70 I = 1,NBIN R(I) = 0.0D0 DR(I) = 0.0D0 SUMM(I) = 0.0 70 CONTINUE C DO 80 I = 1,NREF NSHL = ABS(IBIN(I)) SIGA = SIGMAA(IPS,NSHL) IC = 0 IF (IBIN(I).LT.0) IC = 1 WT = 2.0 - IC X = 2.0*SIGA*EO(I)*EC(I)/ (1.0-SIGA**2) C C ************ CALL SRIN(X,IC,W) C ************ C SUMM(NSHL) = SUMM(NSHL) + W C C---- Omit constants in residual and derivative sums, but C remember to put them in after end of loop. C R(NSHL) = R(NSHL) - DBLE(EO(I)*WT*EC(I)*W) IF (IC.NE.0) THEN C C---- Centric case. (wt=1) C DR(NSHL) = DBLE((EO(I)*EC(I))**2* (W**2-1.0)) + DR(NSHL) ELSE C C---- Non-centric case. C limit as x tends to zero of i1(x)/(x*i0(x)), ie. W/x C for non-centric, is 1/2. avoid potential divide by zero C and substitute if appropriate. include wt=2 C IF (X.LT.1.0E-5) WDX = 0.5 IF (X.GE.1.0E-5) WDX = W/X DR(NSHL) = DBLE((EO(I)*EC(I))**2*4.0* (W**2+WDX-1.0)) + + DR(NSHL) END IF 80 CONTINUE DELMAX = 0.0 OVERM = 0.0 DO 90 I = 1,NBIN IF (NE(I).NE.0) THEN I1 = I + 1 C C---- First add in the constant terms to r and dr C OLD = SIGMAA(IPS,I) R(I) = DBLE(SUMW(I)*OLD) + R(I) DR(I) = DBLE((OLD**2+1.0)/ ((1.0-OLD**2)**2))*DR(I) DR(I) = DBLE(SUMW(I)) + DR(I) C C---- Adjust current estimate of sigmaa; check for convergence C IF (DR(I).NE.0.0) THEN C DEL = REAL(-R(I)/DR(I)) IF (ABS(DEL).GT.DELMAX) DELMAX = ABS(DEL) SIGMAA(IPS,I) = OLD + DEL IF (SIGMAA(IPS,I).GE.1.0) THEN C C---- Sigmaa is restricted to range 0 to 1 C WRITE (ISYSW,FMT=6014) SIGMAA(IPS,I) 6014 FORMAT (' SIGMAA = ',F10.5,' is > 1.0: reset to 0.9', + '99') SIGMAA(IPS,I) = 0.999 END IF W = SUMM(I)/NE(I) OVERM = SUMM(I) + OVERM WRITE (ISYSW,FMT=6016) DLIM(I),DLIM(I1),OLD,DEL, + SIGMAA(IPS,I),W 6016 FORMAT (F7.2,'--',F5.2,4F10.5) ELSE GO TO 140 END IF END IF 90 CONTINUE IF (DELMAX.LT.0.0001) GO TO 110 100 CONTINUE OVERM = OVERM/NREF call ccp4h_summary_beg() WRITE (ISYSW,FMT=6018) OVERM 6018 FORMAT (/,' Refinement of SIGMAA has not Converged -- using ', + 'current values',/' Overall MEAN W is ',F10.5) call ccp4h_summary_end() GO TO 120 110 OVERM = OVERM/NREF call ccp4h_summary_beg() WRITE (ISYSW,FMT=6020) OVERM 6020 FORMAT (/,' Refinement of SIGMAA has CONVERGED',/' Overall ', + 'MEAN FOM is ',F10.5) call ccp4h_summary_end() 120 WRITE (ISYSW,FMT=6022) 6022 FORMAT (' ',/' DATA for SIGMAA PLOT',/' STHOL**2 LN(SIGMAA)') NXY = 0 C DO 130 I = 1,NBIN IF (NE(I).NE.0) THEN c ClemensV Allow for small error.... IF (SIGMAA(IPS,I).GT.0.00001) THEN SIGLOG = LOG(SIGMAA(IPS,I)) WRITE (ISYSW,FMT=6024) ASTL2(I),SIGLOG 6024 FORMAT (2F10.6) NXY = NXY + 1 XY(1,NXY) = ASTL2(I) XY(2,NXY) = SIGLOG END IF END IF 130 CONTINUE WRITE (ISYSW,FMT=6026) 6026 FORMAT (' ln(SigmaA) vs sthol**2 should give straight line with' + ,/' slope = -26.3189*(mean square coordinate error) ', + 'and intercept = 0.5*ln(sigmap/sigman).',/' the firs', + 't and last few points usually deviate and are ignored.', + ///' fit lsq line through points and delete any which ', + 'deviate ',/' by 2*average deviation',/' then refit li', + 'ne. Expected coordinate error given.') C C---- Do sigmaa plot on printer if error estimate needed. C C ****************************************** IF (IERR.EQ.1) CALL PLOTR(XY,NXY,CXY,TITLEX,TITLEY,LEGEND) C ****************************************** C RETURN 140 WRITE (ISYSW,FMT=6028) DLIM(I),DLIM(I1) 6028 FORMAT (' Shift Denominator is 0.0 for ',F6.2,'--',F5.2, + 'A data') C C ********************* CALL CCPERR(1, 'Convergence failure') C ********************* C END IF RETURN C C ***************************************************** 160 CALL CCPERR(1, + ' ***More than ISIZ Reflections Reset Parameter ISIZ **') C ***************************************************** C END C ============================= SUBROUTINE ITRANS(KINT,ICHAR) C ============================= C C---- Translate integers into character strings. C C .. Scalar Arguments .. INTEGER KINT C .. C .. Array Arguments .. CHARACTER ICHAR(5)*1 C .. C .. Local Scalars .. INTEGER I,INT2,J LOGICAL NEG CHARACTER BLANK*1,MINUS*1 C .. C .. Local Arrays .. CHARACTER ITABLE(10)*1 C .. C .. Intrinsic Functions .. INTRINSIC ABS,INT C .. C .. Save statement .. SAVE C .. C .. Data statements .. DATA BLANK/' '/,MINUS/'-'/ DATA ITABLE/'0','1','2','3','4','5','6','7','8','9'/ C .. C NEG = .FALSE. IF (KINT.LT.0) NEG = .TRUE. C C---- Initialize ichar to represent 0, in case it is. C DO 10 I = 1,4 ICHAR(I) = BLANK 10 CONTINUE ICHAR(5) = ITABLE(1) INT2 = ABS(KINT) J = 5 20 CONTINUE C C---- Assume, since it has already been checked, that int C can be printed in 5 characters. C IF (INT2.NE.0) THEN I = INT2 - INT(INT2/10.0)*10 + 1 ICHAR(J) = ITABLE(I) J = J - 1 INT2 = INT(INT2/10.0) GO TO 20 END IF IF (NEG) ICHAR(J) = MINUS END C ==================================================== SUBROUTINE LABEL(VMIN,VMAX,LOWV,INCREM,NTICK,IPOWER) C ==================================================== C C---- From range, determine smallest increment that give C .le. 15 ticks on axis. C C .. Scalar Arguments .. REAL VMAX,VMIN INTEGER INCREM,IPOWER,LOWV,NTICK C .. C .. Local Scalars .. REAL POWER C .. C .. Intrinsic Functions .. INTRINSIC INT,LOG10 C .. C .. Save statement .. SAVE C .. C POWER = LOG10(VMAX-VMIN) C C---- Make sure to correct for fact that ifix of negative C number is not the same as for positive. C IF (POWER.GE.0.0) IPOWER = INT(POWER) - 1 IF (POWER.LT.0.0) IPOWER = INT(POWER) - 2 INCREM = INT((VMAX-VMIN)/ (10.0**IPOWER*15.0)) + 1 C C---- For aesthetic reasons, choose increment of 1,2,5 or 10 C IF (INCREM.GT.2 .AND. INCREM.LT.5) INCREM = 5 IF (INCREM.GT.5 .AND. INCREM.LT.10) INCREM = 10 LOWV = (INT(VMIN/ (10.0**IPOWER*INCREM))+1)*INCREM IF (VMIN.LT.0.0) LOWV = LOWV - INCREM NTICK = INT(VMAX/ (10.0**IPOWER*INCREM)) - LOWV/INCREM + 1 C C---- Check for labels too large or small to print with 5 C characters. ntick of 0 is flag for this. C IF (LOWV.LT.-9999 .OR. (LOWV+NTICK*INCREM).GT.99999) NTICK = 0 END C ==================================================== SUBROUTINE LSQ(X,Y,NOBS,SIG,MWT,A,B,SIGA,SIGB,CHI2,Q) C ==================================================== C C ln(sigmaa) vs sthol**2 should give straight line with C slope = -26.3189*(mean square coord error) and C intercept = 0.5*ln(sigmap/sigman). C the first and last few points usually deviate and are ignored. C to display graph parameter C Calls lsq routine and WRITEs out the numbers associated C with prediced lsq fit C C Author T J Oldfield C C ( York --- January 1988) C C very general lsq program C From "Numerical recipes" C C date January 1988 C C Given a set of nobs points x(i), y(i) with SD sig(i), fit C them to a straight line y=a+bx by minimising x**2. Returned C A,B and their respective probable uncertainties siga and sigb C , the chi-square chi2, and the goodness-of-fit probability C q. IF mwt=0, on input to the routine THEN the sd are assumed C to be unavailable: q is returned as 1.0 and the nobsrmalisation C of chi2 is to unit sd on all the points C C C .. Scalar Arguments .. REAL A,B,CHI2,Q,SIGA,SIGB INTEGER MWT,NOBS C .. C .. Array Arguments .. REAL SIG(NOBS),X(NOBS),Y(NOBS) C .. C .. Local Scalars .. REAL SIGDAT,SS,ST2,SX,SXOSS,SXX,SXY,SY,SYY,T,TMP,WT INTEGER I C .. C .. Intrinsic Functions .. INTRINSIC ABS,SQRT C .. C .. Save statement .. SAVE C .. C SXX = 0 SYY = 0 SXY = 0 SX = 0 SY = 0 ST2 = 0 B = 0 IF (MWT.EQ.0) THEN DO 10 I = 1,NOBS SIG(I) = 1.0 10 CONTINUE END IF SS = 0 DO 20 I = 1,NOBS WT = 1.0/ (SIG(I)**2) SS = SS + WT SXX = (X(I)*WT)**2 + SXX SX = X(I)*WT + SX SXY = (Y(I)*WT)**2* (X(I)*WT) + SXY SYY = (Y(I)*WT)**2 + SYY SY = Y(I)*WT + SY 20 CONTINUE C C---- wt = 1/sig(i)**2 = 1.0 IF sig(i) not set C sx = sum of x*wt C sy = sum of y*wt C ss = sum of wt = nobs IF sig(i) not set C sxx= sum of sq(x*wt) C syy= sum of sq(y*wt) C sxy= sum of x*y*wt**2 C SXOSS = SX/SS C Cdebug write(6,9876)sx,sy,ss,sxx,syy,sxy,mwt C9876 format(' sx sy ss sxx syy sxy mwt',/,6f10.5,i5) C DO 30 I = 1,NOBS T = (X(I)-SXOSS)/SIG(I) ST2 = T*T + ST2 B = Y(I)*T/SIG(I) + B 30 CONTINUE C B = B/ST2 A = (SY-SX*B)/SS SIGA = SQRT((SX*SX/ (SS*ST2)+1.0)/SS) SIGB = SQRT(1.0/ST2) C C---- Q is returned as the correlation coefficient C Q = (SXY- (SX*SY/SS))/SQRT((SXX- (SX*SX/SS))* (SYY- (SY*SY/SS))) CHI2 = 0 C C---- Return Sig(i) as error(y) squared C DO 40 I = 1,NOBS CHI2 = ((Y(I)-A-X(I)*B)/SIG(I))**2 + CHI2 TMP = (Y(I)-A-X(I)*B)/SIG(I) SIG(I) = ABS(TMP) 40 CONTINUE SIGDAT = SQRT(CHI2/ (SS-2.0)) CHI2 = SIGDAT SIGA = SIGA*SIGDAT SIGB = SIGB*SIGDAT Q = SIGDAT C C---- q=gammq(0.5*(nobs-2),0.5*chi2) C C the above is not implemented C END C ============================================================ SUBROUTINE MIRFAZ(IH,NPS,IC,NSHL,FO,A,B,C,D,ALFCOM,WCOM) C ============================================================ C C C .. Scalar Arguments .. REAL A,B,C,D,FO INTEGER IC,NPS,NSHL C .. C .. Array Arguments .. INTEGER IH(3) C .. C .. Arrays in Common .. REAL SMCMIR,SMNMIR INTEGER NCMIR,NNMIR C .. C .. Local Scalars .. REAL ALFCOM,VARM,WCOM INTEGER IOFF C .. C .. External Subroutines .. EXTERNAL MIRFIL C .. C .. Intrinsic Functions .. INTRINSIC COS,SIN C .. C .. Common blocks .. COMMON /STMIR/NCMIR(50),NNMIR(50),SMCMIR(50),SMNMIR(50) C .. C .. Save statement .. SAVE C .. C IF (IC.EQ.0) NNMIR(NSHL) = NNMIR(NSHL) + 1 IF (IC.NE.0) NCMIR(NSHL) = NCMIR(NSHL) + 1 C C---- Fill in mir slots of rdump C C ************************************************ CALL MIRFIL(IH,NPS,IC,NSHL,A,B,C,D,ALFCOM,WCOM,VARM,SMCMIR,SMNMIR) C ************************************************ C C IOFF = NPS*4 + 3 END C =========================================================== SUBROUTINE MIRFIL(IH,NPS,IC,NSHL,A,B,C,D,ALF,W,VARM,SUMMC,SUMMN) C =========================================================== C C C .. Scalar Arguments .. REAL A,B,C,D INTEGER IC,NPS,NSHL C .. C .. Array Arguments .. REAL SUMMC(50),SUMMN(50) INTEGER IH(3) C .. C .. Local Scalars .. REAL ALF,VARINF,VARM,W INTEGER IOFF C .. C .. External Subroutines .. EXTERNAL CENPHA,PHASE C .. C .. Save statement .. SAVE C .. C IOFF = NPS*4 + 3 C IF (IC.NE.0) THEN C C---- Centric C C ***************** CALL CENPHA(A,B,ALF,W) C ***************** C SUMMC(NSHL) = SUMMC(NSHL) + W C RDUMP(IOFF+7) = 0.0 VARM = 0.0 ELSE C C---- Non-centric C C ****************************** CALL PHASE(IH,A,B,C,D,ALF,W,VARINF) C ****************************** C SUMMN(NSHL) = SUMMN(NSHL) + W VARM = VARINF END IF C END C ================================================================= SUBROUTINE PARFAZ(IH,NPS,IC,NSHL,FO,FC,APAR,BPAR,ALFCOM,WCOM, + AOUT, BOUT,AOUTD,BOUTD) C ================================================================= C C C .. Scalar Arguments .. REAL AOUT,AOUTD,BOUT,BOUTD,FO,ALFCOM,WCOM INTEGER IC,NPS,NSHL C .. C .. Array Arguments .. INTEGER IH(3) REAL APAR(3),BPAR(3),FC(3),ALFPAR(3),WPAR(3),VARP(3) C .. C .. Arrays in Common .. REAL CDPAR,SMCCP,SMCPAR,SMNCP,SMNPAR,VARCP,VARPAR,WTPAR, + SIGMAA,SIGMAN,SIGMAP,ACHANG INTEGER NCPAR,NNPAR C .. C .. Local Scalars .. REAL ATOT,BTOT,DLUZ,VARINF,WTDEN C .. C .. External Subroutines .. EXTERNAL CENPHA,PARFIL,PHASE,UNBIAS C .. C .. Intrinsic Functions .. INTRINSIC COS,SIN C .. C .. Common blocks .. COMMON /SIGS/SIGMAA(3,50),SIGMAN(3,50),SIGMAP(3,50),ACHANG(3) COMMON /STPAR/NCPAR(50),NNPAR(50),SMCPAR(3,50),SMNPAR(3,50), + VARPAR(3,50),CDPAR(3,50),SMCCP(50),SMNCP(50),VARCP(50), + WTPAR(3,50) C .. C .. Save statement .. SAVE C .. C IF (IC.EQ.0) NNPAR(NSHL) = NNPAR(NSHL) + 1 IF (IC.NE.0) NCPAR(NSHL) = NCPAR(NSHL) + 1 ATOT = 0.0 BTOT = 0.0 WTDEN = 0.0 C C---- Fill in partial structure slots of rdump. C C ******************************************************** CALL PARFIL(IH,NPS,IC,NSHL,APAR,BPAR,ALFPAR,WPAR,VARP, + ATOT,BTOT,WTDEN,SMCPAR,SMNPAR,VARPAR) c CALL PARFIL(IH,NPS,IC,NSHL,FC,APAR,BPAR,RDUMP,ATOT,BTOT,WTDEN, c + SMCPAR,SMNPAR,VARPAR) C ******************************************************** C C---- Fill in rest of rdump, and calculate map coeffs. C C---- Avoid some work if only one partial structure C IF (NPS.NE.1) THEN C C **************************** IF (IC.NE.0) CALL CENPHA(ATOT,BTOT,ALFCOM,WCOM) C **************************** C IF (IC.EQ.0) THEN C C **************************************** CALL PHASE(IH,ATOT,BTOT,0.0,0.0,ALFCOM,WCOM,VARINF) C **************************************** C C RDUMP(3) = VARINF VARCP(NSHL) = VARCP(NSHL) + VARINF END IF ELSE C ALF = RDUMP(5) C W = RDUMP(6) C IF (IC.EQ.0) RDUMP(3) = RDUMP(7) ALFCOM = ALFPAR(1) WCOM = WPAR(1) END IF C IF (IC.EQ.0) SMNCP(NSHL) = SMNCP(NSHL) + WCOM IF (IC.NE.0) SMCCP(NSHL) = SMCCP(NSHL) + WCOM C RDUMP(1) = ALF C RDUMP(2) = W AOUT = WCOM*FO*COS(ALFCOM) BOUT = WCOM*FO*SIN(ALFCOM) C C write(6,'(a,/,3i3,10f10.4))') C 1 ' in parfaz1', C 1 ih,apar(1),bpar(1),alfpAR(1),wpAR(1),varp(1), C 1 wtden,alfcom,wcom,aout,bout C---- Remove partial structure bias from map C C -- comments from Randy.. C NOW WE HAVE COMBINED PHASE WEIGHTED FO. THIS IS ALL WE NEED C FOR CENTRIC OPTION 3 COEFFICIENTS. FOR NON-CENTRIC OPTION 3, C REMOVE PARTIAL STRUCTURE BIAS WITH UNBIAS. FOR NON-CENTRIC C OPTION 4, UNBIAS JUST KEEPS SOME STATISTICS. FOR OPTION 4, C TAKE THE DIFFERENCE WITH PARTIAL STRUCTURE 1 HERE. C DLUZ = SIGMAA(1,NSHL) * SQRT(SIGMAN(1,NSHL)/SIGMAP(1,NSHL)) AOUTD = AOUT - DLUZ*FC(1)*COS(ALFPAR(1)) BOUTD = BOUT - DLUZ*FC(1)*SIN(ALFPAR(1)) C C **************************************************** c IF (IC.EQ.0) CALL UNBIAS(ALF,RDUMP,WTDEN,NPS,NSHL,AOUT,BOUT, c + WTPAR,CDPAR) CALL UNBIAS(IC,ALFCOM,FC,ALFPAR,VARP,WTDEN,NPS,NSHL, + AOUT,BOUT,WTPAR,CDPAR) C **************************************************** C write(6,'(a,/,3i3,10f10.4))') C 1 ' in parfaz2',ih,apar(1),bpar(1),alfpar(1),wpar(1),varp(1), C 1 wtden,alfcom,wcom,aout,bout RETURN C END C ============================================================== SUBROUTINE PARFIL(IH,NPS,IC,NSHL,APAR,BPAR,ALFPAR,WPAR,VARP, + ATOT,BTOT, WTDEN,SUMMC,SUMMN,VARPAR) C ============================================================== C C---- Fill in rdump for partial structure. C C .. Scalar Arguments .. REAL ATOT,BTOT,WTDEN INTEGER IC,NPS,NSHL C .. C .. Array Arguments .. REAL APAR(3),BPAR(3),ALFPAR(3),WPAR(3),VARP(3), + SUMMC(3,50),SUMMN(3,50),VARPAR(3,50) INTEGER IH(3) C .. C .. Local Scalars .. INTEGER I C .. C .. External Subroutines .. EXTERNAL CENPHA,PHASE C .. C .. Save statement .. SAVE C .. C DO 10 I = 1,NPS ATOT = APAR(I) + ATOT BTOT = BPAR(I) + BTOT IF (IC.NE.0) THEN C C ***************************** c CALL CENPHA(APAR(I),BPAR(I),ALF,W) CALL CENPHA(APAR(I),BPAR(I),ALFPAR(I),WPAR(I)) C ***************************** C c RDUMP(IOFF+4) = 0.0 SUMMC(I,NSHL) = SUMMC(I,NSHL) + WPAR(I) ELSE C C ********************************************** CALL PHASE(IH,APAR(I),BPAR(I),0.0,0.0,ALFPAR(I),WPAR(I), + VARP(I)) C ********************************************** C c RDUMP(IOFF+4) = VARINF WTDEN = WTDEN + VARP(I) SUMMN(I,NSHL) = SUMMN(I,NSHL) + WPAR(I) VARPAR(I,NSHL) = VARPAR(I,NSHL) + VARP(I) END IF C c RDUMP(IOFF+2) = ALF c RDUMP(IOFF+3) = W 10 CONTINUE C write(6,'(a,3i3,7f10.3))') C 1 ' in parfil',ih,apar(1),bpar(1),alfpAR(1),wpAR(1),varp(1),alf,w END C ====================================================== SUBROUTINE PHASE(IH,AHL,BHL,CHL,DHL,ALPHB,FIGM,VARINF) C ====================================================== C C ALPHB derived from HL coefficients: ALPHB in radians C C .. Scalar Arguments .. REAL AHL,ALPHB,BHL,CHL,DHL,FIGM,VARINF C .. C .. Array Arguments .. INTEGER IH(3) C .. C .. Scalars in Common .. REAL ARGMAX INTEGER NSTEP C .. C .. Arrays in Common .. REAL COSA,DLIM,ETAB,SINA C .. DOUBLE PRECISION EXPARD C .. Local Scalars .. REAL ARG,ARGHI,COSAA,EXPARG,NHL,PROB,SINAA,TWOPI,WIDTH INTEGER I,IND,ICNT LOGICAL GOOD C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN2,INT,LOG,REAL,SQRT C .. C .. Common blocks .. COMMON /TABS/NSTEP,COSA(120),SINA(120),ETAB(4001),ARGMAX,DLIM(51) C .. C .. Save statement .. SAVE C .. C .. Data statements .. DATA TWOPI/6.283185/ C .. C C---- Use hendrickson-lattman coefficients to determine non-centric C phase and figure of merit, in addition to variation of C information. C C---- Width is 2*pi/nstep, ie. width in radians of rectangle C used for numerical integration (rectangular rule) C WIDTH = TWOPI/REAL(NSTEP) ICNT = 0 NHL = 0.0 10 CONTINUE GOOD = .TRUE. ARGHI = 0.0 VARINF = 0.0 COSAA = 0.0 SINAA = 0.0 PROB = 0.0 DO 20 I = 1,NSTEP ARG = COSA(I)*AHL + NHL + SINA(I)*BHL + + (1.0-SINA(I)**2*2.0)*CHL + DHL*2.0*SINA(I)*COSA(I) C C---- exp(arg) is unnormalized phase probability density C (hendrickson and lattman). C etab will rarely be exceeded in practice. we skip angles C with large negative arguments -- these are insignificant C for the sums, and approximating by anything other than 0 C leads to problems with values summed to get varinf. if C there are any values of arg for which etab does not go C high enough, store the highest value and repeat with a C value of nhl that prevents overshooting of etab. negative C arguments are accomodated by taking the inverse. C IF (ARG.GE.-ARGMAX) THEN IF (ARG.GT.ARGMAX) THEN GOOD = .FALSE. IF (ARG.GT.ARGHI) ARGHI = ARG END IF C C---- Avoid some work if we won't be calculating phases etc. C IF (GOOD) THEN IND = INT(ABS(ARG)*50.0+1.5) IF(IND.GT.4001) IND = 4001 IF(IND.LT.1) IND = 1 IF (ARG.GE.0.0) EXPARD = DBLE(ETAB(IND)*WIDTH) IF (ARG.LT.0.0) EXPARD = DBLE(WIDTH/ETAB(IND)) EXPARG = 0.0 IF(EXPARD.GE.1.D-8) EXPARG = SNGL(EXPARD) PROB = PROB + EXPARG COSAA = COSA(I)*EXPARG + COSAA SINAA = SINA(I)*EXPARG + SINAA VARINF = EXPARG*ARG + VARINF END IF END IF 20 CONTINUE IF (.NOT.GOOD) THEN C C---- Phase probability curve has been clipped, so adjust nhl C to make arghi less than the maximum allowed for etab. C this is only a pseudo-normalization, but it suffices: C strictly speaking, nhl should be log of factor required C to normalize integrated probability density to 1. C ICNT = ICNT +1 NHL = ARGMAX - 0.01 - ARGHI C Avoid infinite loop - set FIGM = 0 if this looks like happening IF(ICNT.LT.10)GO TO 10 SINAA = 0.0 COSAA = 0.0 PROB = 1.0 END IF C C---- Prob is integral of probability density, which will not C in general integrate to 1 because we are picking nhl to C avoid jumping out of etab, not to make the integral work C out. C FIGM = SQRT((SINAA/PROB)**2+ (COSAA/PROB)**2) IF (FIGM.NE.0.0) ALPHB = ATAN2(SINAA,COSAA) IF (FIGM.EQ.0.0) ALPHB = 0.0 IF (ALPHB.LT.0.0) ALPHB = ALPHB + TWOPI C Cdebug write(6,9876)ih,ahl,bhl,chl,dhl,alp,figm C9876 format(' in phase h k l a b c d alphb figm',3i3,4f8.3,5x,2f8.3) C C---- Variation of information (varinf) is a measure of the amount C of information added to a posterior distribution compared to C a prior probability density. with an uninformative prior, it C is a measure of the information content of the probability C density. varinf is given by: C integral(0 to 2pi) [p(alpha)*ln(p(alpha)/p0(alpha))*dalpha] C ("information theory with applications", silviu guiasu (1977)) C in this case, p0(alpha) = 1.0/2pi, ie. an uninformative prior. C VARINF = MAX(0.0,(LOG(TWOPI/PROB) + VARINF/PROB)) END C ================================================= SUBROUTINE PLOTR(XY,NXY,CXY,TITLEX,TITLEY,LEGEND) C ================================================= C C---- Given an array of x,y pairs, this s/r fills in an array which, C when printed in "a" format, constitutes a 2-d graph of C the points, properly scaled and nicely labelled. the C arguments supplied are, in order, the array of x,y pairs, C number of points, the character to represent the points, C title(in 20a1 format) to be printed under the x axis, title C to be printed along the y axis, and legend to be printed under C the x axis title. C C---- Variables for lsq line fit - test which points should C be excluded to give a good slope used to estimate C coordinate error C C .. Scalar Arguments .. INTEGER NXY CHARACTER CXY*1 C .. C .. Array Arguments .. REAL XY(2,50) CHARACTER LEGEND(20)*1,TITLEX(20)*1,TITLEY(20)*1 C .. C .. Scalars in Common .. REAL CONV,DELST2,SMAX,SMAXCB,SMIN,SMINCB,ST2MIN,WMAX,WMAX2 INTEGER IERR,IOPT,ISYSR,ISYSW,NBIN,NCOL,NMON,NPS,MIR2,IRFSIG C .. C .. Local Scalars .. REAL A,B,CHI2,CORRN,Q,RMSERR,SIGA,SIGB,SIGRMS,XMAX,XMIN,XRANGE, + YMAX,YMIN,YRANGE INTEGER I,I1,I2,ICYC,IFLAG,INCX,INCY,IPOWRX,IPOWRY,ISCALE,J,K, + KINT,LOWX,LOWY,MWT,NCYC,NOBS,NOBS0,NTICKX,NTICKY CHARACTER BAR*1,BLANK*1,CORNER*1,DASH*1,LE*1,ONE*1,STAR*1,TICK*1 C .. C .. Local Arrays .. REAL SIG(50),XLSQ(50),YLSQ(50) CHARACTER ICHAR(5)*1,PLOT(129,60)*1 C .. C .. External Subroutines .. EXTERNAL CCPERR,ITRANS,LABEL,LSQ C .. C .. Intrinsic Functions .. INTRINSIC ABS,INT,SQRT C .. C .. Common blocks .. COMMON /CNTROL/NBIN,ST2MIN,DELST2,NPS,MIR2,IOPT,SMAX,IERR,WMAX, + WMAX2,SMIN,SMINCB,SMAXCB,IRFSIG COMMON /IO/ISYSR,ISYSW,NCOL,CONV,NMON C .. C .. Save statement .. SAVE C .. C .. Data statements .. C C---- Choose most suitable characters available on printer. C DATA BLANK/' '/,BAR/'|'/,DASH/'-'/,ONE/'1'/,CORNER/'+'/ DATA TICK/'+'/,LE/'E'/,STAR/'*'/ C .. C C---- SLOPE = -26.3189*(MEAN SQUARE COORD ERROR) C CORRN = -1.0/26.3189 NCYC = 10 ICYC = 1 DO 10 I = 1,NXY XLSQ(I) = XY(1,I) YLSQ(I) = XY(2,I) SIG(I) = 1.0 10 CONTINUE 20 CONTINUE C IF (ICYC.EQ.1) THEN NOBS = NXY IF (NOBS.GT.50) GO TO 160 END IF WRITE (6,FMT='(''1'')') J = 0 DO 30 I = 1,NOBS IF (SIG(I).NE.0.0) THEN J = J + 1 XLSQ(J) = XLSQ(I) YLSQ(J) = YLSQ(I) SIG(J) = 1.0 END IF 30 CONTINUE NOBS = J cejd Don't try to estimate coordinate error when you are wanging! IF (IERR.EQ.1) THEN CEJD IF (NOBS.NE.NOBS0) THEN NOBS0 = NOBS MWT = 1 C C ************************************************ CALL LSQ(XLSQ,YLSQ,NOBS,SIG,MWT,A,B,SIGA,SIGB,CHI2,Q) C ************************************************ C C---- SLOPE = -26.3189*(MEAN SQUARE COORD ERROR) AND C RMSERR = SQRT(ABS(B*CORRN)) SIGRMS = -SIGB*CORRN WRITE (6,FMT=*) ' ' WRITE (6,FMT=*) ' ' IF (B.LE.0.0) call ccp4h_summary_beg() WRITE (6,FMT=*) ' ' WRITE (6,FMT='(A)') ' From plot of ln(SigmaA) vs sthol**2 :' WRITE (6,FMT=*) ' ' WRITE (6,FMT='(A,F15.6)') ' Least square line gradient = ' + ,B WRITE (6,FMT='(A,F15.6)') ' Uncertainty in gradient = ' + ,SIGB WRITE (6,FMT='(A,F15.6)') ' x intercept = ' + ,A WRITE (6,FMT='(A,F15.6)') ' Uncertainty in intercept = ' + ,SIGA WRITE (6,FMT=*) ' ' C C---- Problems in estimating coordinate error if b is positive C IF (B.GT.0.00001) THEN WRITE (6,FMT='(A,/,A)') ' ** PROBLEM!! B IS POSITIVE **', + ' *** NOT POSSIBLE TO ESTIMATE COORDINATE ERROR .***' DO I = 1,NOBS IF (SIG(I).GT.2.0*CHI2) SIG(I) = 0.0 ENDDO ELSE RMSERR = SQRT(B*CORRN) SIGRMS = -SIGB*CORRN WRITE (6,FMT='(A,F15.6)') ' Rms coordinate error = ',RMSERR WRITE (6,FMT='(A,F15.6)') + ' UNCERTAINTY IN RMS ERROR SQUARED= ',SIGRMS call ccp4h_summary_end() C C---- Returned Sig(i) = error(ylsq(i) ** 2) C call ccp4h_graph_beg(0,0) WRITE (6,FMT=6038) 6038 FORMAT(/, $ ' $TABLE: Estimation of rms coordinate error:',/, $ ' $GRAPHS',/, $ ' :ln(SigmaA) vs sthol**2 ', $ ' :N:2,3:',/, $ ' $$ Bin sthol**2 ln(SigmaA) error iflag $$', $ /,' $$',/) DO 40 I = 1,NOBS IFLAG = 0 IF (SIG(I).GT.2.0*CHI2) IFLAG = 1 WRITE (6,FMT=6000) I,XLSQ(I),YLSQ(I),SIG(I),IFLAG 6000 FORMAT (I8,3F10.4,I4) IF (IFLAG.EQ.1) SIG(I) = 0.0 40 CONTINUE WRITE(6,FMT='(/,A,//)') ' $$' call ccp4h_graph_end() END IF ICYC = ICYC + 1 IF (ICYC.LT.NCYC) GO TO 20 END IF END IF XMIN = 1.0E30 XMAX = -1.0E30 YMIN = 1.0E30 YMAX = -1.0E30 C C---- Determine outer bounds of x and y values C IF (NXY.GE.1) THEN DO 50 I = 1,NXY IF (XY(1,I).LT.XMIN) XMIN = XY(1,I) IF (XY(1,I).GT.XMAX) XMAX = XY(1,I) IF (XY(2,I).LT.YMIN) YMIN = XY(2,I) IF (XY(2,I).GT.YMAX) YMAX = XY(2,I) 50 CONTINUE XRANGE = XMAX - XMIN YRANGE = YMAX - YMIN IF (XRANGE.NE.0.0 .AND. YRANGE.NE.0.0) THEN C C---- Add 5% margin on each edge C XMAX = 0.05*XRANGE + XMAX XMIN = XMIN - 0.05*XRANGE YMAX = 0.05*YRANGE + YMAX YMIN = YMIN - 0.05*YRANGE XRANGE = 1.1*XRANGE YRANGE = 1.1*YRANGE C C---- Determine how axes are to be labelled. C C **************************************** CALL LABEL(XMIN,XMAX,LOWX,INCX,NTICKX,IPOWRX) CALL LABEL(YMIN,YMAX,LOWY,INCY,NTICKY,IPOWRY) C **************************************** C C---- Now initialize blanks and axes for character array C representing plot of graph C DO 70 J = 1,60 DO 60 I = 1,129 PLOT(I,J) = BLANK IF (I.EQ.9 .AND. J.LT.55) PLOT(I,J) = BAR IF (J.EQ.55 .AND. I.GT.9) PLOT(I,J) = DASH 60 CONTINUE 70 CONTINUE PLOT(1,1) = ONE PLOT(9,55) = CORNER DO 80 I = 1,20 I1 = I + 40 I2 = I + 20 PLOT(I1,58) = TITLEX(I) PLOT(I1,60) = LEGEND(I) PLOT(2,I2) = TITLEY(I) 80 CONTINUE C C---- If x axis labels can fit into 5 characters, label x axis. C IF (NTICKX.NE.0) THEN KINT = LOWX - INCX DO 100 I = 1,NTICKX KINT = KINT + INCX J = INT((10.0**IPOWRX*KINT-XMIN)/XRANGE*120.0) + 9 IF (J.GE.9 .AND. J.LE.125) THEN PLOT(J,55) = TICK C C ****************** CALL ITRANS(KINT,ICHAR) C ****************** C DO 90 I2 = 1,5 IF (ICHAR(I2).NE.BLANK) THEN PLOT(J,56) = ICHAR(I2) J = J + 1 END IF 90 CONTINUE END IF 100 CONTINUE C C---- If labels have been multiplied by a factor of 10, indicate C this factor on the x axis title. C IF (IPOWRX.NE.0) THEN ISCALE = -IPOWRX C C ******************** CALL ITRANS(ISCALE,ICHAR) C ******************** C PLOT(62,58) = STAR PLOT(63,58) = ONE PLOT(65,58) = LE I = 66 C DO 110 J = 1,5 IF (ICHAR(J).NE.BLANK) THEN PLOT(I,58) = ICHAR(J) I = I + 1 END IF 110 CONTINUE END IF END IF C C---- Now label y axis C IF (NTICKY.NE.0) THEN KINT = LOWY - INCY DO 130 I = 1,NTICKY KINT = KINT + INCY J = INT((YMAX-10.0**IPOWRY*KINT)/YRANGE*55.0) IF (J.GT.0 .AND. J.LE.55) THEN PLOT(9,J) = TICK C C ****************** CALL ITRANS(KINT,ICHAR) C ****************** C DO 120 I2 = 1,5 K = I2 + 3 PLOT(K,J) = ICHAR(I2) 120 CONTINUE END IF 130 CONTINUE C C---- Indicate multiple of 10. C IF (IPOWRY.NE.0) THEN ISCALE = -IPOWRY C C ******************** CALL ITRANS(ISCALE,ICHAR) C ******************** C PLOT(2,42) = STAR PLOT(2,43) = ONE PLOT(2,45) = LE I = 46 C DO 140 J = 1,5 IF (ICHAR(J).NE.BLANK) THEN PLOT(2,I) = ICHAR(J) I = I + 1 END IF 140 CONTINUE END IF END IF C C---- Finally ready to map points into array. C DO 150 I = 1,NXY J = INT((XY(1,I)-XMIN)/XRANGE*120.0) + 9 K = INT((YMAX-XY(2,I))/YRANGE*55.0) PLOT(J,K) = CXY 150 CONTINUE C C---- Plot is ready for printer. C WRITE (ISYSW,FMT=6004) PLOT 6004 FORMAT (59 (129A1,/),129A1) END IF END IF C RETURN C C ************************************************* 160 CALL CCPERR(1, + 'Too many observations -- reset NDATA in code') C ************************************************* C END C ================ SUBROUTINE SFOUT C ================ C C---- Write out map coefficients for maps using phase C information from only one partial structure (iopt = 1 or 2). C accumulate an excessive amount of statistics. C C C .. Parameters .. INTEGER MCOLS PARAMETER (MCOLS=500) C .. C .. Scalars in Common .. REAL ARGMAX,CONV,DELST2,SMAX,ST2MIN,WMAX,WMAX2 INTEGER IAC1,IAC2,IAC3,IAO,IAO2,ICOA,ICOB,ICOC,ICOD,IERR,IFC1, + IFC2,IFC3,IFO,IFS,IOPT,ISYSR,ISYSW,IW,IW2,NBIN,NCOL, + NMON,NPS,NSTEP,MIR2,IRFSIG C .. C .. Arrays in Common .. REAL COSA,DLIM,ETAB,SIGMAA,SIGMAN,SIGMAP,SINA,ACHANG C .. C .. Local Scalars .. REAL ALF,AVEMCT,AVEMNT,EC,EO,EPSI,FC,FO,FOUT,OVERM,SSQ, + STOL2,VARINF,W,X INTEGER I,I1,IAC,IC,ICHK,IFC,INDX,ISYSAB,J,NCTOT,NNTOT,NREF,NSHL LOGICAL EOF C .. C .. Local Arrays .. REAL ADATA(MCOLS),AVEMC(50),AVEMN(50),R(11) INTEGER HIST(10,50),HISTT(10),IH(3),NC(50),NN(50) LOGICAL LOGMSS(MCOLS) C .. C .. External Subroutines .. EXTERNAL CENTR,EPSLON,LRREFL,LWREFL,SRIN C .. C .. Intrinsic Functions .. INTRINSIC INT,MOD,NINT,REAL,SQRT C .. C .. Common blocks .. COMMON /CNTROL/NBIN,ST2MIN,DELST2,NPS,MIR2,IOPT,SMAX,IERR,WMAX, + WMAX2,SMIN,SMINCB,SMAXCB,IRFSIG COMMON /IO/ISYSR,ISYSW,NCOL,CONV,NMON COMMON /POSN/IFO,IFS,IAO,IW,ICOA,ICOB,ICOC,ICOD,IFC1,IAC1,IFC2, + IAC2,IFC3,IAC3,IAO2,IW2,ICOA2,ICOB2,ICOC2,ICOD2 COMMON /SIGS/SIGMAA(3,50),SIGMAN(3,50),SIGMAP(3,50),ACHANG(3) COMMON /TABS/NSTEP,COSA(120),SINA(120),ETAB(4001),ARGMAX,DLIM(51) C .. C .. Save statement .. SAVE C .. C VARINF = 0.0 C C---- Make header for DUMP file, if it was requested C CONV = 3.14159/180.0 C DO 20 I = 1,NBIN NC(I) = 0 AVEMC(I) = 0.0 NN(I) = 0 AVEMN(I) = 0.0 DO 10 J = 1,10 HIST(J,I) = 0 10 CONTINUE 20 CONTINUE C C---- Ejd 13-Sep-1988 I think that Fo and Fc are assumed C to be on the same scale which is not C necessarily true for our MTZ file input C - Grrrrrr C C---- Calculate fourier coefficients. for 2fo-fc type map, C use 2*m*fobs - d*fcalc, alpha calc for non-centric, m*fobs C for centric. this removes model bias. since m*fobs, C alpha calc is best estimate for fobs, alpha true, use C m*fobs - fcalc, alpha calc for difference map C C IFO=lookup(5) C IFS=lookup(6) C IAO=lookup(7) C IW=lookup(8) C ICOE=lookup(9) C IFC1=lookup(13) C IAC1=lookup(14) C IFC2=lookup(15) C IAC2=lookup(16) C IFC3=lookup(17) C IAC3=lookup(18) C C see above C SMAX = 4.0*DELST2*NBIN + 4.0*ST2MIN C C---- Here IPS must equal 0 C IFC = IFC1 IAC = IAC1 NREF = 0 C WRITE (ISYSW,FMT=6000) 6000 FORMAT (///) WRITE (ISYSW,FMT=*) + ' SFOUT -ref no - h k l Eo Ec Nshl Ic Epsi' 30 CONTINUE C C---- This reads one piece of data from the MTZ file C and returns it in ADATA - (500 REAL array) C C *********************** CALL LRREFL(1,SSQ,ADATA,EOF) C *********************** C IF (EOF) GO TO 40 C CALL LRREFM(1,LOGMSS) C C---- Set 6 data items after the last input item = the Magic value C CALL EQUAL_MAGIC(1,ADATA(NCOL+1),6) C C---- Check sd(f) and resolution limits C IF (SSQ.GT.SMAX .OR. SSQ.LT.SMIN) GO TO 30 IF (LOGMSS(IFO) .AND. LOGMSS(IFC)) GO TO 35 IF (.NOT.LOGMSS(IFS)) THEN IF (ADATA(IFS).LT.0.000001) GO TO 30 ENDIF C NREF = NREF + 1 IH(1) = NINT(ADATA(1)) IH(2) = NINT(ADATA(2)) IH(3) = NINT(ADATA(3)) C C---- Prepare for output C FO = 0.0 FC = 0.0 ALF = 0.0 IF (.NOT.LOGMSS(IFO)) FO = ADATA(IFO) IF (.NOT.LOGMSS(IFC)) FC = ADATA(IFC) IF (.NOT.LOGMSS(IAC)) THEN ALF = ADATA(IAC)*CONV ELSE GOTO 35 ENDIF STOL2 = SSQ/4.0 NSHL = INT((STOL2-ST2MIN)/DELST2) + 1 IF (NSHL.LT.1) NSHL = 1 IF (NSHL.GT.NBIN) NSHL = NBIN C C ********************** CALL EPSLON(IH,EPSI,ISYSAB) CALL CENTR(IH,IC) C ********************** C EO = FO/SQRT(SIGMAN(1,NSHL)*EPSI) EC = FC/SQRT(SIGMAP(1,NSHL)*EPSI) X = SIGMAA(1,NSHL)*2.0*EO*EC/ (1.0-SIGMAA(1,NSHL)**2) cv C. Vonrhein Jul 5 1999 c c If we don't have FP we still want to output coefficients c IF (LOGMSS(IFO)) X = 1. cv C C ************ C Big X W tends to 1.0; X=0 W=0 CALL SRIN(X,IC,W) C ************ C Modify W if desired W = W*ACHANG(1) DLUZ = SIGMAA(1,NSHL) * SQRT(SIGMAN(1,NSHL)/SIGMAP(1,NSHL)) C ICHK = MOD(NREF,NMON) IF (ICHK.EQ.1) WRITE (ISYSW,FMT=6002) NREF,IH,EO,EC,NSHL,IC, + EPSI 6002 FORMAT (7X,I7,3X,3I3,4X,2F6.3,2I4,F4.0) C C---- Fill in histogram of figures of merit for each shell. C merge centric and non-centric data for this. C INDX = INT(W*10.0) + 1 IF (INDX.GT.10) INDX = 10 HIST(INDX,NSHL) = HIST(INDX,NSHL) + 1 IF (IC.EQ.1) THEN C C---- Here for centric data C C CENTRIC DATA: EITHER M*FO OR M*FO-D*FC C FOUT = W*FO cv C. Vonrhein Jul 5 1999 c c if we don't have FP set it to -D*Fc c IF (LOGMSS(IFO)) FOUT = -DLUZ*FC c FOUTD = FOUT -DLUZ*FC NC(NSHL) = NC(NSHL) + 1 AVEMC(NSHL) = AVEMC(NSHL) + W ELSE C C---- Here to get estimate of fobs, alpha true for non-centric C have to remove model bias for 2fo-fc type coeffs, but not C for difference map coeffs. C C NON-CENTRIC DATA: EITHER 2*M*FO-D*FC OR M*FO-D*FC C C FOUT = (2.0*W*EO-SIGMAA(1,NSHL)*EC)*SQRT(SIGMAN(1,NSHL)*EPSI) C DFOUT = W*FO FOUT = 2.0*W*FO - DLUZ*FC FOUTD = W*FO - DLUZ*FC NN(NSHL) = NN(NSHL) + 1 AVEMN(NSHL) = AVEMN(NSHL) + W END IF C C---- All data again. is it a difference map? C C---- Ejd 13-Sep-1988 I think that Fo and Fc are assumed C to be on the same scale which is not necessarily C true for our MTZ file input - Grrrrrr C C DFOUT = DFOUT - SQRT(SIGMAN(1,NSHL)*EPSI)*EC C C---- LRREFL output C DATA NLPRGO,LSPRGO/10,'MFODC','PHMFODC','2MFODFC','PH2MFODFC', C + 'PHCMB','WCMB', 'HLAC','HLBC','HLCC','HLDC', C C ADATA(NCOL+1) = W C ADATA(NCOL+2) = DFOUT C ADATA(NCOL+3) = FOUT C ADATA(NCOL+4) = ALF/CONV c ADATA(NCOL+1) = ABS(FOUTD) c ADATA(NCOL+2) = ALF/CONV c IF(FOUTD.LT.0.0) ADATA(NCOL+2) = ALF/CONV + 180.0 c ADATA(NCOL+3) = ABS(FOUT) c ADATA(NCOL+4) = ALF/CONV c IF(FOUT.LT.0.0) ADATA(NCOL+4) = ALF/CONV + 180.0 c ADATA(NCOL+5) = W c ADATA(NCOL+6) = ALF/CONV IF (.NOT.LOGMSS(IFC)) THEN IF (.NOT. (LOGMSS(IFO))) ADATA(NCOL+1) = FOUTD ADATA(NCOL+2) = FOUT IF (LOGMSS(IFO)) ADATA(NCOL+2) = -FOUT C = DLUZ*FC for acentrics, and 0.0 for centrics END IF C ADATA(NCOL+3) = W cv C. Vonrhein Jul 5 1999 c c make sure that the weight is set to zero in cases with no FP c IF (LOGMSS(IFO)) ADATA(NCOL+3) = 0. cv C C *************** 35 CALL LWREFL(1,ADATA) C *************** C GO TO 30 40 NNTOT = 0 AVEMNT = 0.0 NCTOT = 0 AVEMCT = 0.0 WRITE (ISYSW,FMT=6004) 6004 FORMAT (' Figure of Merit Statistics for Partial Structure',/' ', + ' Centric Data Non-Centric ALL Data', + /' D Limits N N N ') C DO 50 I = 1,NBIN IF ((NC(I)+NN(I)).NE.0) THEN NREF = NC(I) + NN(I) OVERM = (AVEMN(I)+AVEMC(I))/NREF NNTOT = NN(I) + NNTOT AVEMNT = AVEMN(I) + AVEMNT IF (NN(I).NE.0) AVEMN(I) = AVEMN(I)/NN(I) NCTOT = NC(I) + NCTOT AVEMCT = AVEMC(I) + AVEMCT IF (NC(I).NE.0) AVEMC(I) = AVEMC(I)/NC(I) I1 = I + 1 WRITE (ISYSW,FMT=6006) DLIM(I),DLIM(I1),NC(I),AVEMC(I),NN(I), + AVEMN(I),NREF,OVERM 6006 FORMAT (F7.2,'--',F5.2,3 (I7,F8.4)) END IF C 50 CONTINUE NREF = NCTOT + NNTOT IF (NREF.NE.0) OVERM = (AVEMCT+AVEMNT)/NREF IF (NCTOT.NE.0) AVEMCT = AVEMCT/NCTOT IF (NNTOT.NE.0) AVEMNT = AVEMNT/NNTOT WRITE (ISYSW,FMT=6008) NCTOT,AVEMCT,NNTOT,AVEMNT,NREF,OVERM 6008 FORMAT (' Overall ',3 (I7,F8.4)) C DO 60 I = 1,11 R(I) = REAL(I-1)/10.0 60 CONTINUE C DO 70 I = 1,10 HISTT(I) = 0 70 CONTINUE WRITE (ISYSW,FMT=6010) R 6010 FORMAT (' ',/' Figure of Merit Histograms (all DATA)',/' D Li', + 'mits ',F4.1,10F5.1) C DO 90 I = 1,NBIN DO 80 J = 1,10 HISTT(J) = HISTT(J) + HIST(J,I) 80 CONTINUE C I1 = I + 1 WRITE (ISYSW,FMT=6012) DLIM(I),DLIM(I1), (HIST(J,I),J=1,10) 6012 FORMAT (F7.2,'--',F5.2,10I5) 90 CONTINUE C WRITE (ISYSW,FMT=6014) HISTT 6014 FORMAT (' Overall ',10I5) call ccp4h_summary_beg() WRITE(6,395) 395 FORMAT(/,' Weighted Fourier coefficients written to HKLOUT.',/) call ccp4h_summary_end() RETURN C C---- DUMP everything of relevance, if requested. put 0.0 for C variation of information to indicate that this was not C calculated. same output as in s/r DUMPer, except it stops C at fc. C END C ==================== REAL FUNCTION SIM(X) C ==================== C C---- Calculate sim and srinivasan non-centric figure of merit C as i1(x)/i0(x), where i1 and i0 are the modified 1st and zero C order bessel functions. C references: sim, g. a. (1960) acta cryst. 13, 511-512; C srinivasan, r. (1966) acta cryst. 20, 143-144; C abramowitz & stegun, handbook of mathematical functions, 378. C this routine was obtained from w. kabsch. C C .. Scalar Arguments .. REAL X C .. C .. Local Scalars .. REAL T C .. C .. Intrinsic Functions .. INTRINSIC ABS,SIGN C .. C .. Save statement .. SAVE C .. C T = ABS(X)/3.75 C IF (T.GT.1.0) THEN T = 1.0/T SIM = ((((((((0.01787654-T*0.00420059)*T+ (-0.02895312))*T+ + 0.02282967)*T+ (-0.01031555))*T+0.00163801)*T+ + (-0.00362018))*T+ (-0.03988024))*T+0.39894228)* + SIGN(1.0,X)/ ((((((((-0.01647633+T*0.00392377)*T+ + 0.02635537)*T+ (-0.02057706))*T+0.00916281)*T+ + (-0.00157565))*T+0.00225319)*T+0.01328592)*T+0.39894228) ELSE T = T*T SIM = ((((((T*0.00032411+0.00301532)*T+0.02658733)*T+ + 0.15084934)*T+0.51498869)*T+0.87890594)*T+0.5)*X/ + ((((((T*0.0045813+0.0360768)*T+0.2659732)*T+1.2067492)*T+ + 3.0899424)*T+3.5156229)*T+1.0) END IF END C ========================== SUBROUTINE SRIN(X,IC,WCAL) C ========================== C C---- Calculate figure of merit, given x. formula differs for C centric and non-centric. srinivasan (acta cryst 20:143(1966)) C differs from woolfson (acta cryst 9:804(1956)) and sim (acta C cryst 13:511(1960)) only in the expression for x, which in C the present case includes contributions from missing structure C and model error. C for non-centric reflections, use the function sim which is C the ratio of the modified 1st and zero order bessel functions. C C .. Scalar Arguments .. REAL WCAL,X INTEGER IC C .. C .. External Functions .. REAL SIM EXTERNAL SIM C .. C .. Intrinsic Functions .. INTRINSIC TANH C .. C .. Save statement .. SAVE C .. C IF (IC.EQ.0) WCAL = SIM(X) IF (IC.NE.0) WCAL = TANH(X/2.0) C END C =============================== SUBROUTINE STATC(NPS,NBIN,DLIM) C =============================== C C---- Put out stats for combined phases. C C C .. Scalar Arguments .. INTEGER NBIN,NPS C .. C .. Array Arguments .. REAL DLIM(51) C .. C .. Scalars in Common .. REAL CONV INTEGER ISYSR,ISYSW,NCOL,NMON C .. C .. Arrays in Common .. REAL CDMIR,CDPC,SMMCOM,SMCMC,SMCPC,SMNCOM,SMNMC,SMNPC,VARCOM, + VARMIR,VARPC,WEIGHT INTEGER NCCOM,NNCOM C .. C .. Local Scalars .. INTEGER I,I1,J,NCTOT,NNTOT C .. C .. Common blocks .. COMMON /IO/ISYSR,ISYSW,NCOL,CONV,NMON COMMON /STCOM/NCCOM(50),NNCOM(50),SMCPC(3,50),SMNPC(3,50), + VARPC(3,50),CDPC(3,50),SMCMC(50),SMNMC(50),VARMIR(50), + CDMIR(50),SMMCOM(50),SMNCOM(50),VARCOM(50),WEIGHT(3,50) C .. C .. Save statement .. SAVE C .. C IF (NPS.NE.0) THEN NCTOT = 0 NNTOT = 0 C DO 10 I = 1,NBIN NCTOT = NCCOM(I) + NCTOT NNTOT = NNCOM(I) + NNTOT 10 CONTINUE C IF ((NCTOT+NNTOT).NE.0) THEN WRITE (ISYSW,FMT=6000) 6000 FORMAT (' ',/' Reflections with Combined Phase Information') IF (NCTOT.NE.0) THEN WRITE (ISYSW,FMT=6002) (I,I=1,NPS) 6002 FORMAT (' Centric DATA Mean Figures of MERIT', + /41X,'Partial Structures',/' D Limits N C', + 'ombined MIR',3I10) C DO 30 I = 1,NBIN IF (NCCOM(I).NE.0) THEN SMMCOM(I) = SMMCOM(I)/NCCOM(I) SMCMC(I) = SMCMC(I)/NCCOM(I) DO 20 J = 1,NPS SMCPC(J,I) = SMCPC(J,I)/NCCOM(I) 20 CONTINUE I1 = I + 1 WRITE (ISYSW,FMT=6004) DLIM(I),DLIM(I1),NCCOM(I), + SMMCOM(I),SMCMC(I), (SMCPC(J,I),J=1,NPS) 6004 FORMAT (F7.2,'--',F5.2,I6,5F10.4) END IF 30 CONTINUE END IF IF (NNTOT.NE.0) THEN WRITE (ISYSW,FMT=6006) (I,I=1,NPS) 6006 FORMAT (' ',/' Non-Centric DATA',/' ', + ' Mean Figures of MERIT',/41X,'Partial Struct', + 'ures',/' D Limits N Combined MIR ', + 3I10) C DO 50 I = 1,NBIN IF (NNCOM(I).NE.0) THEN SMNCOM(I) = SMNCOM(I)/NNCOM(I) DO 40 J = 1,NPS SMNPC(J,I) = SMNPC(J,I)/NNCOM(I) 40 CONTINUE SMNMC(I) = SMNMC(I)/NNCOM(I) I1 = I + 1 WRITE (ISYSW,FMT=6004) DLIM(I),DLIM(I1),NNCOM(I), + SMNCOM(I),SMNMC(I), (SMNPC(J,I),J=1,NPS) END IF 50 CONTINUE WRITE (ISYSW,FMT=6008) (I,I=1,NPS) 6008 FORMAT (' Mean Variation of Information',/ + ' Partial Structures' + ,/' D Limits Combined MIR ',3I10) C DO 70 I = 1,NBIN IF (NNCOM(I).NE.0) THEN VARCOM(I) = VARCOM(I)/NNCOM(I) DO 60 J = 1,NPS VARPC(J,I) = VARPC(J,I)/NNCOM(I) 60 CONTINUE VARMIR(I) = VARMIR(I)/NNCOM(I) I1 = I + 1 WRITE (ISYSW,FMT=6010) DLIM(I),DLIM(I1),VARCOM(I), + VARMIR(I), (VARPC(J,I),J=1,NPS) 6010 FORMAT (F7.2,'--',F5.2,5F10.4) END IF 70 CONTINUE C WRITE (ISYSW,FMT=6012) (I,I=1,NPS) 6012 FORMAT (' Mean Cosines of Differences from Combine', + 'd PHASE',/' Partial Struct', + 'ures',/' D Limits MIR',3I10) C DO 90 I = 1,NBIN IF (NNCOM(I).NE.0) THEN DO 80 J = 1,NPS CDPC(J,I) = CDPC(J,I)/NNCOM(I) 80 CONTINUE CDMIR(I) = CDMIR(I)/NNCOM(I) I1 = I + 1 WRITE (ISYSW,FMT=6010) DLIM(I),DLIM(I1),CDMIR(I), + (CDPC(J,I),J=1,NPS) END IF 90 CONTINUE WRITE (ISYSW,FMT=6014) (I,I=1,NPS) 6014 FORMAT (' Mean Partial Structure Weights', + /' Partial Structures',/' D Limi', + 'ts ',I8,2I10) C DO 110 I = 1,NBIN IF (NNCOM(I).NE.0) THEN DO 100 J = 1,NPS WEIGHT(J,I) = WEIGHT(J,I)/NNCOM(I) 100 CONTINUE I1 = I + 1 WRITE (ISYSW,FMT=6010) DLIM(I),DLIM(I1), + (WEIGHT(J,I),J=1,NPS) END IF 110 CONTINUE END IF END IF END IF C END C =========================== SUBROUTINE STATM(NBIN,DLIM) C =========================== C C---- Put out stats for mir phase info only. C C C .. Scalar Arguments .. INTEGER NBIN C .. C .. Array Arguments .. REAL DLIM(51) C .. C .. Scalars in Common .. REAL CONV INTEGER ISYSR,ISYSW,NCOL,NMON C .. C .. Arrays in Common .. REAL SMCMIR,SMNMIR INTEGER NCMIR,NNMIR C .. C .. Local Scalars .. INTEGER I,I1,ITEST C .. C .. Common blocks .. COMMON /IO/ISYSR,ISYSW,NCOL,CONV,NMON COMMON /STMIR/NCMIR(50),NNMIR(50),SMCMIR(50),SMNMIR(50) C .. C .. Save statement .. SAVE C .. C ITEST = 0 C DO 10 I = 1,NBIN ITEST = NCMIR(I) + ITEST + NNMIR(I) 10 CONTINUE C IF (ITEST.NE.0) THEN WRITE (ISYSW,FMT=6000) 6000 FORMAT (' ',/' Reflections with only MIR PHASES',/'0 ', + ' Centric DATA Non-Centric',/' D Limits ', + ' N N ') C DO 20 I = 1,NBIN IF ((NCMIR(I)+NNMIR(I)).NE.0) THEN IF (NCMIR(I).NE.0) SMCMIR(I) = SMCMIR(I)/NCMIR(I) IF (NNMIR(I).NE.0) SMNMIR(I) = SMNMIR(I)/NNMIR(I) I1 = I + 1 WRITE (ISYSW,FMT=6002) DLIM(I),DLIM(I1),NCMIR(I),SMCMIR(I), + NNMIR(I),SMNMIR(I) 6002 FORMAT (F7.2,'--',F5.2,2 (I7,F8.4)) END IF 20 CONTINUE C END IF C END C =============================== SUBROUTINE STATP(NPS,NBIN,DLIM) C =============================== C C---- Put out stats for partial structure phase info only. C C C .. Scalar Arguments .. INTEGER NBIN,NPS C .. C .. Array Arguments .. REAL DLIM(51) C .. C .. Scalars in Common .. REAL CONV INTEGER ISYSR,ISYSW,NCOL,NMON C .. C .. Arrays in Common .. REAL CDPAR,SMCCP,SMCPAR,SMNCP,SMNPAR,VARCP,VARPAR,WTPAR INTEGER NCPAR,NNPAR C .. C .. Local Scalars .. INTEGER I,I1,J,NCTOT,NNTOT C .. C .. Common blocks .. COMMON /IO/ISYSR,ISYSW,NCOL,CONV,NMON COMMON /STPAR/NCPAR(50),NNPAR(50),SMCPAR(3,50),SMNPAR(3,50), + VARPAR(3,50),CDPAR(3,50),SMCCP(50),SMNCP(50),VARCP(50), + WTPAR(3,50) C .. C .. Save statement .. SAVE C .. C IF (NPS.NE.0) THEN NCTOT = 0 NNTOT = 0 C DO 10 I = 1,NBIN NCTOT = NCPAR(I) + NCTOT NNTOT = NNPAR(I) + NNTOT 10 CONTINUE IF ((NCTOT+NNTOT).NE.0) THEN WRITE (ISYSW,FMT=6000) 6000 FORMAT (' Reflections with only Partial Structure PHASES') C IF (NCTOT.NE.0) THEN C WRITE (ISYSW,FMT=6002) (I,I=1,NPS) 6002 FORMAT (' Centric DATA Mean Figures of MERIT',/ + ' Partial Structures', + /' D Limits N Combined ',I6,2I8) C DO 30 I = 1,NBIN IF (NCPAR(I).NE.0) THEN SMCCP(I) = SMCCP(I)/NCPAR(I) C DO 20 J = 1,NPS SMCPAR(J,I) = SMCPAR(J,I)/NCPAR(I) 20 CONTINUE C I1 = I + 1 WRITE (ISYSW,FMT=6004) DLIM(I),DLIM(I1),NCPAR(I), + SMCCP(I), (SMCPAR(J,I),J=1,NPS) 6004 FORMAT (F7.2,'--',F5.2,I7,4F8.4) END IF 30 CONTINUE END IF C IF (NNTOT.NE.0) THEN WRITE (ISYSW,FMT=6006) (I,I=1,NPS) 6006 FORMAT (' Non-Centric Mean Figures of MERIT',/ + ' Partial Structures', + /' D Limits N Combined ',I6,2I8) C DO 50 I = 1,NBIN IF (NNPAR(I).NE.0) THEN SMNCP(I) = SMNCP(I)/NNPAR(I) DO 40 J = 1,NPS SMNPAR(J,I) = SMNPAR(J,I)/NNPAR(I) 40 CONTINUE I1 = I + 1 WRITE (ISYSW,FMT=6004) DLIM(I),DLIM(I1),NNPAR(I), + SMNCP(I), (SMNPAR(J,I),J=1,NPS) END IF 50 CONTINUE IF (NPS.NE.1) THEN WRITE (ISYSW,FMT=6008) (I,I=1,NPS) 6008 FORMAT (' Mean Variation of Information', + /' Partial Structures', + /' D Limits Combined',I8,2I10) C DO 70 I = 1,NBIN IF (NNPAR(I).NE.0) THEN VARCP(I) = VARCP(I)/NNPAR(I) DO 60 J = 1,NPS VARPAR(J,I) = VARPAR(J,I)/NNPAR(I) 60 CONTINUE I1 = I + 1 WRITE (ISYSW,FMT=6010) DLIM(I),DLIM(I1),VARCP(I), + (VARPAR(J,I),J=1,NPS) 6010 FORMAT (F7.2,'--',F5.2,4F10.4) END IF 70 CONTINUE C WRITE (ISYSW,FMT=6012) (I,I=1,NPS) 6012 FORMAT (' Mean Cosines of Differences from Combine', + 'd Phase',/' Partial Structures', + /' D Limits ',I9,2I10) C DO 90 I = 1,NBIN IF (NNPAR(I).NE.0) THEN DO 80 J = 1,NPS CDPAR(J,I) = CDPAR(J,I)/NNPAR(I) 80 CONTINUE I1 = I + 1 WRITE (ISYSW,FMT=6010) DLIM(I),DLIM(I1), + (CDPAR(J,I),J=1,NPS) END IF 90 CONTINUE C WRITE (ISYSW,FMT=6014) (I,I=1,NPS) 6014 FORMAT (' Mean Partial Structure Weights', + /' Partial Structures',/' D Li', + 'mits ',I9,2I10) C DO 110 I = 1,NBIN IF (NNPAR(I).NE.0) THEN DO 100 J = 1,NPS WTPAR(J,I) = WTPAR(J,I)/NNPAR(I) 100 CONTINUE I1 = I + 1 WRITE (ISYSW,FMT=6010) DLIM(I),DLIM(I1), + (WTPAR(J,I),J=1,NPS) END IF 110 CONTINUE END IF END IF END IF END IF C END C ====================================== SUBROUTINE STINIT(NPS,NBIN,HIST,HISTT) C ====================================== C C---- Just to make things a little easier to read, do all the C variable initialization for combined phase stats here. C C C .. Scalar Arguments .. INTEGER NBIN,NPS C .. C .. Array Arguments .. INTEGER HIST(10,50),HISTT(10) C .. C .. Arrays in Common .. REAL CDMIR,CDPAR,CDPC,SMMCOM,SMCCP,SMCMC,SMCMIR,SMCPAR,SMCPC, + SMNCOM,SMNCP, + SMNMC,SMNMIR,SMNPAR,SMNPC,VARCOM,VARCP,VARMIR,VARPAR,VARPC, + WEIGHT,WTPAR INTEGER NCCOM,NCMIR,NCPAR,NNCOM,NNMIR,NNPAR C .. C .. Local Scalars .. INTEGER I,J C .. C .. Common blocks .. COMMON /STCOM/NCCOM(50),NNCOM(50),SMCPC(3,50),SMNPC(3,50), + VARPC(3,50),CDPC(3,50),SMCMC(50),SMNMC(50),VARMIR(50), + CDMIR(50),SMMCOM(50),SMNCOM(50),VARCOM(50),WEIGHT(3,50) COMMON /STMIR/NCMIR(50),NNMIR(50),SMCMIR(50),SMNMIR(50) COMMON /STPAR/NCPAR(50),NNPAR(50),SMCPAR(3,50),SMNPAR(3,50), + VARPAR(3,50),CDPAR(3,50),SMCCP(50),SMNCP(50),VARCP(50), + WTPAR(3,50) C .. C .. Save statement .. SAVE C .. C DO 20 I = 1,NBIN NCCOM(I) = 0 NNCOM(I) = 0 SMCMC(I) = 0.0 SMNMC(I) = 0.0 VARMIR(I) = 0.0 CDMIR(I) = 0.0 SMMCOM(I) = 0.0 SMNCOM(I) = 0.0 VARCOM(I) = 0.0 NCPAR(I) = 0 NNPAR(I) = 0 SMCCP(I) = 0.0 SMNCP(I) = 0.0 VARCP(I) = 0.0 NCMIR(I) = 0 NNMIR(I) = 0 SMCMIR(I) = 0.0 SMNMIR(I) = 0.0 IF (NPS.NE.0) THEN DO 10 J = 1,NPS SMCPC(J,I) = 0.0 SMNPC(J,I) = 0.0 VARPC(J,I) = 0.0 CDPC(J,I) = 0.0 WEIGHT(J,I) = 0.0 SMCPAR(J,I) = 0.0 SMNPAR(J,I) = 0.0 VARPAR(J,I) = 0.0 CDPAR(J,I) = 0.0 WTPAR(J,I) = 0.0 10 CONTINUE END IF 20 CONTINUE DO 40 I = 1,10 HISTT(I) = 0 DO 30 J = 1,NBIN HIST(I,J) = 0 30 CONTINUE 40 CONTINUE END C ================= SUBROUTINE TABLES C ================= C C .. Scalars in Common .. REAL ARGMAX INTEGER NSTEP C .. C .. Arrays in Common .. REAL COSA,DLIM,ETAB,SINA C .. C .. Local Scalars .. REAL ANG,ARG,D2R,STEP INTEGER I C .. C .. Intrinsic Functions .. INTRINSIC COS,EXP,REAL,SIN C .. C .. Common blocks .. COMMON /TABS/NSTEP,COSA(120),SINA(120),ETAB(4001),ARGMAX,DLIM(51) C .. C .. Save statement .. SAVE C .. C .. Data statements .. DATA D2R/1.745329E-02/ C .. C C---- Set up lookup tables for cos, sin and exp. C STEP = 360.0/NSTEP DO 10 I = 1,NSTEP ANG = (I-1)*D2R*STEP COSA(I) = COS(ANG) SINA(I) = SIN(ANG) 10 CONTINUE C C---- Set up table of exp(arg) for arg from 0 to 80. inverses will be C used for negative arguments. step size of 0.02 ensures that C maximum error is 1 percent, which is quite sufficient. since C exp(80)=5.5e34, there will be no overflow even on vaxes and C the range of possible values will far exceed what is necessary C for the phase probability curve integration. C DO 20 I = 1,4001 ARG = REAL(I-1)/50.0 ETAB(I) = EXP(ARG) 20 CONTINUE ARGMAX = 80.0 END C =============================================================== SUBROUTINE UNBIAS(IC,ALFCOM,FC,AC,VARC,WTDEN,NPS,NSHL,AOUT,BOUT, + WEIGHT, COSDEL) C =============================================================== C Phases all in radians C C---- Calculate map coefficients to remove model bias from C combined phase map. C for mir phase info, m*|fo|, alpha mir ~ fo C for partial structure, m*|fo|, alpha par ~ fo/2 + d*fc/2 C assume that combined phase will give a value that varies C between these as a linear function of w, which varies from C 0 to 1 with amount of partial structure influence. C we can include several partial structures to get: C m*|fo|, alpha comb ~ [(2-sum(wI))/2 *fo + sum (wI*dI/2)*fcI] C then solve for fo (remember these are vectors with phases): C fo ~ (m*|fo|,alpha com - sum[(wI*dI/2)*fcI]) / (1-sum(wI)/2) C for the weight, we use the variation of information for the C phase probability from that source, divided by the sum of C the varinf's from all of the sources (so the w's add up to 1.) C C .. Scalar Arguments .. REAL ALFCOM,AOUT,BOUT,WTDEN INTEGER IC,NPS,NSHL C .. C .. Array Arguments .. C REAL FC(3),AC(3),WC(3),COSDEL(3,50),RDUMP(22),WEIGHT(3,50) REAL FC(3),AC(3),VARC(3),COSDEL(3,50),WEIGHT(3,50) C .. C .. Arrays in Common .. REAL SIGMAA,SIGMAN,SIGMAP,ACHANG C .. C .. Local Scalars .. REAL DENOM,DLUZ,SUMWT INTEGER I, IFO C .. C .. Intrinsic Functions .. INTRINSIC COS,SIN,SQRT C .. C .. Common blocks .. COMMON /SIGS/SIGMAA(3,50),SIGMAN(3,50),SIGMAP(3,50),ACHANG(3) C .. C .. Save statement .. SAVE C .. C SUMWT = 0.0 C cv C. Vonrhein Jul 5 1999 c use this subroutine also for centrics with no Fobs cv IF ((ABS(AOUT).GT.0.).OR.(ABS(BOUT).GT.0.)) THEN IFO = 1 ELSE IFO = 0 END IF C IF ((IFO.EQ.1).AND.(IC.NE.0)) RETURN C DO 10 I = 1,NPS c IOFF = (I-1)*4 + 3 c FC = RDUMP(IOFF+1) c ALF = RDUMP(IOFF+2) DLUZ = SQRT(SIGMAN(I,NSHL)/SIGMAP(I,NSHL))*SIGMAA(I,NSHL) c WT = RDUMP(IOFF+4) IF (WTDEN.NE.0.0) VARC(I) = VARC(I)/WTDEN SUMWT = SUMWT + VARC(I) C Write(6,'(a,i3,12f10.4)') C 1 ' in unbias1',i,dluz,fc(i),ac(i),VARc(i),wtden,sumwt,aout,bout cv IF (IC.EQ.0) THEN AOUT = AOUT - VARC(I)/2.0*DLUZ*FC(I)*COS(AC(I)) BOUT = BOUT - VARC(I)/2.0*DLUZ*FC(I)*SIN(AC(I)) ELSE AOUT = AOUT - DLUZ*FC(I)*COS(AC(I)) BOUT = BOUT - DLUZ*FC(I)*SIN(AC(I)) END IF c Write(6,'(a,i3,12f10.4)') c 1 ' in unbias2',i,dluz,fc(i),ac(i),VARc(i),wtden,sumwt,aout,bout WEIGHT(I,NSHL) = WEIGHT(I,NSHL) + VARC(I) COSDEL(I,NSHL) = COS(AC(I)-ALFCOM) + COSDEL(I,NSHL) 10 CONTINUE IF (IC.EQ.0) THEN DENOM = 1.0 - SUMWT/2.0 ELSE DENOM = 1.0 END IF AOUT = AOUT/DENOM BOUT = BOUT/DENOM C cv C. Vonrhein Jul 5 1999 c invert sign for cases without Fobs cv IF (IFO.EQ.0) THEN AOUT = -AOUT BOUT = -BOUT END IF C Write(6,'(a,i3,12f10.4)') C 1 ' in unbias3',i,dluz,fc(1),ac(1),VARc(1),wtden,sumwt,aout,bout END