C C This code is distributed under the terms and conditions of the C CCP4 licence agreement as `Part ii)' software. See the conditions C in the CCP4 manual for a copyright statement. C C************************************************************* C C Modified on Mon Jul 14 09:02:44 BST 1997 C for harvest sub calls for deposition to PDB C C Program Converted to MTZ from LCF data file format C C Date: Mon Jan 13 16:37:31 GMT 1992 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 PROGRAM TRUNCATE C ================= C CTitle:TRUNCATE REFLECTION DATA - TRUNCATE C****************************************** C C CDate:Tue, 11 Dec 90 16:27 GMT CFrom:Eleanor Dodson C C-- TRUNCATE TRUNCATE.FOR 7/3/86 JWC C C C=====MTZ CONVERSION C *************** C C TRUNCATE C C C---- DESCRIPTION C =========== C C The program reads a file of averaged intensities (SCALA output) C and writes a file containing mean amplitues and anomalous differences C with or without truncating the data. C C Truncates intensities (assumed normally distributed) to amplitudes C Calculates mean and anomalous difference C C C K.S.WILSON & S.FRENCH OXFORD 1976/7 C C THIS VERSION PHILIP EVANS APRIL 1977 C C C Modified 14-November-1990 C Uses new Symmetry options. No need for CENT card - Symmetry read as C nspgrp or SPACE GROUP NAME. C C C Modified to standard PARSER input KEYWORDING 1-2-1989 C see new document in [PUBLIC.XTAL.DOC.UTX]TRUNCATE.UTX C C modified for new AGROVATA output, A.C.Bloomer Jun 1981 C further modified for triclinic reflections JMS1 Sep 198 C LCF version for VAX January 1982 PRE C C STOP program if any reflection is outside the 60th bin, since this can C lead to reflections being incorrectly rejected. 31/10/86 A.G.W.L. C C C C Replaced linear interpolation with one based on cubic spline. C (This makes significant difference when interpolating SIGMA C i.e. the average intensity for given resolution. In general the C distribution of observed intensities gets closer to theoretical C values. (In my test cases, an average discrepancy of ca 2% of reflections C went down to <1%.) On the other hand, the effect of the spline C on Wilson distribution is very small. The distribution seems to be C sampled finely enough in the tables, so it makes little difference C how you interpolate.) C C If necessary, adjust the width of bin to cover the resolution limits C evenly. This applies for the case when the width of bin was set C explicitely on input and there is a tail-end of a few reflections C in the last bin. C C Increase default number of bins from 30 to 60. C C If .lt. MINREF (40) reflections in any one bin then: C if RANGES was not set explicitely on input then: C - increase the width of bins and try again. C if RANGES was set: C - print a warning C C Stop if width of bin greater than WDTHMX (0.03) C C Some tidying up. Eg, calculate width of bin between resolution limits C and not between resmax and infinity. C C WRR October 1992 C C------------------------------------------------------------------------ C C Input MTZ file is the ouput file from AGROVATA with 1 record C per reflection having format C C H K L Imean SD IW+ SD I+ SD IW- SD I- SD C C where the primed values are used for the anomalous differences C Centric data has only 5 columns set C acentric up to 13 columns set C C Output MTZ file has format C H K L Fmean SD DeltaF SD ISYM C C where ISYM=0 normally but for those HKLs where all data C refer to I+ or to I- ISYM is set to 1 or 2 indicating tha C FMEAN value is biased. C After optional truncation (and/or SQRT) of as many of C Imean I+ I- I+' I-' as are present , giving the correspond C FMEAN=0.5*(F+ + F-) for acentric C FMEAN= F+ or F- if one & only one C FMEAN=Fmean for centric C DELTAF= (F+' - F-') C C C This file includes subroutines and functions: C S SETCLL(C) C F STHLQ4(INHKL) C S NEGIAS(XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ) C S NEGICS(XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ) C S SCNDER(X,Y,N,YP1,YPN,Y2) calculate 2nd derivatives for cubic spline C S SPLINT(XA,YA,Y2A,N,X,Y) returns interpolated values C C C EXTERNAL SUBROUTINES & FUNCTIONS ARE USED FOR: C DISTINGUISHING CENTRIC AND ACENTRIC REFLECTIONS C CENTR EPSABS IN SYMLIB C C================================================================== C C March 1998: FALLOFF program from Yorgo Modis incorporated C by Martyn Winn. Calculates the falloff of mean F-values in C 3 orthogonal directions, and produces plot file illustrating C this. C C------------------------------------------------------------------------ C C PROGRAM_FUNCTION C C The program TRUNCATE reads a reflection data file C of averaged intensities (SCALA output) and outputs an MTZ C reflection data file containing F and DeltaFanom values. The C input intensities are assumed to follow a normal distribution C with the standard deviations, i.e. negative observations must C have been preserved. The truncation procedure used was devised C by French and Wilson (reference 1) and is based on Bayesian C statistics. The F's are calculated using the prior knowledge of C Wilson's distributions for acentric or centric data (calculated C in shells of reciprocal space in a first pass through the data) C and the mean intensity and standard deviation values. The F's C output are all positive and follow Wilson's distribution. The C truncation procedure has little effect on reflections larger C than 3 standard deviations but should give significantly better C values for the weak data than those obtained by merely taking the C square root of the intensities and setting negative intensities C to zero. Reflections of less than minus four standard deviations C are rejected. C C The following warnings should be heeded: C C a) Do not truncate data more than once, e.g. do not merge C truncated data with untruncated data and then truncate again. C C b) The standard deviations are crucial to the procedure C and the standard deviation analysis from SCALA should be checked. C C c) The procedure should not be used on data which has been forced C to be positive e.g. from ordinate analysis measurements on a C diffractometer or where negative observations have been set to zero. C C REFERENCES C C 1) French G.S. and Wilson K.S. Acta. Cryst. (1978), A34 , 517. C C C IMPLICIT NONE C C .. Parameters .. INTEGER NLPRGI PARAMETER (NLPRGI=16) INTEGER ncols PARAMETER (ncols=19) INTEGER MAXSETS PARAMETER (MAXSETS=ncols) INTEGER MAXHIS PARAMETER (MAXHIS=30) INTEGER MAXSER PARAMETER (MAXSER=60) INTEGER MAXSYM PARAMETER (MAXSYM=96) INTEGER NPARM PARAMETER (NPARM=200) C..maximum number of bins INTEGER KNBINS,KNBNS2,KNBS20 PARAMETER (KNBINS=60) PARAMETER (KNBNS2=2*KNBINS, KNBS20=20*KNBINS) C..maximum allowed width of bin REAL WDTHMX PARAMETER (WDTHMX=0.03) C..minimum allowed number of reflection in any one bin INTEGER MINREF PARAMETER (MINREF=40) C .. C .. Scalars in Common .. REAL RANGE INTEGER ISOUT,KANOM,LABFLG,LBOTKN,LENTIT,LUNIN,LUNOUT,MTITOT, + MTZBPR,MTZERR,MTZFLG,MTZIN,MTZOUT,MTZPRT,NBATX, + NHSOUT,NLAUE,NREFLX,NSPGRP,NSYM,NSYMP, + NUMCOL LOGICAL CSYMM,KYCELL,KYHIST,KYSYMM,KYTITL,KYRESO,LBOLOG, $ LNOTRN,MTZEOF,RNGSET CHARACTER LATTYP*1,DATEMT*8,LAUNAM*10,NAMSPG*10,PGNAME*10, + VERSNX*10,TITIN*70,TITOUT*70,LBOLIN*400 C .. C .. Arrays in Common .. REAL CELL,CELMTZ,CELOUT,RNGMTZ,RSYM,RSYMT INTEGER LBOBEG,LBOEND,MTZLOK CHARACTER CTPRGO*1,CTPRGI*1,LSPRGO*30, + LSPRGI*30,HISOUT*80 C .. C .. Local Scalars .. REAL DMAX,F,FMAX,FMIN,FN,RMAXI,S,SDF,SDJ,SIGMA,SIGMAX, + SIGMIN,SMAX,SMIN,SQWT,SS,UPPER,WEIGHT,XF,XJ INTEGER I,IAPPND,IBATCH,IC,ICEN,IFAIL,IOUT,IREJ,IS, + ISYMB,ISYSAB,J,JKD,JLOOP,JLOP,K,LDUM, + NCIP,NCSP,NFRAC,NRANGE,NREF,NLPRGO, + NREJ,NT,JJ,JDO,JDO220,JDO240,JDO250,JDO260,JDO270, $ JDO280,JDO290,JDO300,JDO310,JDO320,MINBIN REAL TOTIN,TOTII,TOTIS,TOTIR,TOTFI,TOTFS,TOTFR, $ AIMEAN,SIGIMEAN,AFMEAN,SIGFMEAN,ZERO LOGICAL LTRAP,LOGOUT C .. C .. Local Arrays .. REAL FA(10),FF(KNBINS),FFA(19),FLOWER(KNBINS),FRACS(KNBINS,10,2), + FRACT(10,2),RECIN(16),RECOUT(19),SD(KNBINS), + SDIA(19),SSA(19),THA(10),THC(10),XIA(19) REAL FFD2(KNBINS),RNGCEN(KNBINS) INTEGER INHKL(3),N(KNBINS),NPOSS(KNBINS),NTOT(2),NW(KNBINS,2) REAL SUMIIS(KNBINS),SUMISS(KNBINS),SUMIRS(KNBINS),SUMINS(KNBINS), $ SUMFIS(KNBINS),SUMFSS(KNBINS),SUMFRS(KNBINS) LOGICAL LOGMSS(NCOLS),QISNAN C C Systematic absences INTEGER MAXSYSABS PARAMETER (MAXSYSABS=500) INTEGER IHKLSYSABS(3,MAXSYSABS),NSYSABS REAL FSIGSYSABS(2,MAXSYSABS) C C FALLOFF integer nmin,nmax real aa,bb,cc,alpha,beta,gamma,rad,fvol real ast,bst,cst,alphast,betast,gammast real transf(3,3),sthmax,sthmin,hc(3),ang(3) real somdir(3,KNBINS),somov(KNBINS),somsddir(3,KNBINS), + somsdov(KNBINS),sthstap,cone integer numdir(3,KNBINS),numov(KNBINS) logical KYFALL,PLOTX,PLOTFALL COMMON /COORD/nmin,nmax,sthstap,somdir,somov,somsddir,somsdov, + numdir,numov COMMON /FALL/KYFALL,PLOTX,PLOTFALL,cone C .. C .. External Functions .. INTEGER IDOT EXTERNAL IDOT C .. C .. External Subroutines .. EXTERNAL ASUSET,CCPDAT,CCPERR,CCPFYP,CCPDPN,CCPRCS,CENTR,CENTRIC, + EPSLN,EPSLON,HKLLIM,KEYIN,LKYOUT,LRASSN,LRCELL,LRCLOS,LRINFO, + LROPEN,LRREFF,LRREFM,LRREWD,LRSYMI,LRSYMM,LRTITL,LWASSN, + LWCELL,LWCLOS,LWHIST,LWOPEN,LWREFL,LWSYMM,LWTITL,NEGIAS, + NEGICS,SETCLL,SCNDER,SPLINT,LWHSTL,QISNAN,CCP4H_INIT, + CCP4H_TOC_BEG,CCP4H_TOC_ENT,CCP4H_TOC_END,CCP4H_PRE_BEG, + CCP4H_PRE_END,CCP4H_RULE,CCP4H_LINK,CCP4H_LINK_KEY, + CCP4H_HEADER,CCP4H_GRAPH_BEG,CCP4H_GRAPH_END, + CCP4H_SUMMARY_BEG,CCP4H_SUMMARY_END C .. C .. Intrinsic Functions .. INTRINSIC MAX,MIN,NINT,REAL,SQRT C .. C .. Common blocks .. C COMMON /AGRSYM/NSPGRP,NSYM,NSYMP,NLAUE,CSYMM, + RSYM(4,4,MAXSYM),RSYMT(4,4,MAXSYM) COMMON /CHRMTZ/LSPRGI(16),CTPRGI(16), + LSPRGO(ncols),CTPRGO(ncols),HISOUT(MAXHIS), + VERSNX,TITIN,TITOUT,DATEMT,LBOLIN COMMON /CHRSYM/PGNAME,NAMSPG,LAUNAM,LATTYP COMMON /INOUT/LUNIN,LUNOUT,ISOUT COMMON /KEYISC/KANOM,LNOTRN COMMON /KEYRAR/CELL(6) COMMON /KEYRSC/RANGE,NRANGE COMMON /WRKMTZ/MTZIN,MTZOUT,MTZERR,MTZBPR,MTZPRT,MTZFLG,MTITOT, + NBATX,NHSOUT,NUMCOL,NREFLX,LENTIT,LABFLG,MTZEOF, + KYTITL,KYHIST,KYCELL,KYSYMM,KYRESO,MTZLOK(NLPRGI), + CELMTZ(6),CELOUT(6),RNGMTZ(2,16), + LBOBEG(NPARM),LBOEND(NPARM),LBOLOG,LBOTKN COMMON /RNG/RNGSET C .. COMMON /APLSCL/ SCALEF, LAPSCL REAL SCALEF LOGICAL LAPSCL C .. REAL SI,SN,SR,SW,XY,XYAV,FFS,PI,DELDSS,F000, + VOL,RR,CELLAS,RVOL,FRAC INTEGER NRFUNQ,NATMS,NELEC,ISMAX COMMON /PREWIL/ SI(100),SN(100),SR(100),SW(100), + NRFUNQ(100),XY(2,100),XYAV(2,100), + FFS(500),PI,DELDSS,NATMS,NELEC,F000, + ISMAX,VOL,RR(3,3,6),CELLAS(6),RVOL,FRAC C C---- Wilson additions here C INTEGER NBIN,NMON,IPLFSQ,NK REAL RMAXWL,RMINWL,SMINWL,SMAXWL,SMINSC,SMAXSC,RMINSC,RMAXSC, $ VPA COMMON /WILKEY/ NBIN,NMON,RMAXWL,RMINWL,SMINWL,SMAXWL,SMINSC, + SMAXSC,RMINSC,RMAXSC,VPA, + IPLFSQ,NK LOGICAL SCBLOG,KYTWIN INTEGER IWITCH,NREFWL REAL WILSC(2),WILBT(2) COMMON /LOGWIL/ IWITCH,SCBLOG,KYTWIN C CHARACTER NWW*2 CHARACTER NWW*4 COMMON /WILCHR/ NWW(20) INTEGER NC,KatmWghts REAL FormFactAcoeff,FormFactBcoeff COMMON /WILFRM/ NC(20),KatmWghts(20),FormFactAcoeff(5,20), + FormFactBcoeff(4,20) C C----- C C Storage for /**M for M = 1, MAXMMT C SMMEMM = sum(E**M) for M, resolution, acentric/centric C Take E as sqrt(I/) First and third moment,/**0.5 I**1.5>/**1.5 C For theoretical distribution. C when k is integer ( moments of I) kth-moment of I = k!. C when k is 0.5, 1.5, 2.5 etc (moments of E E3 E5..) C moment = 0.5** sqrt(PI), 1.5*0.5*sqrt(PI), 2.5*1.5*0.5*sqrt(PI) etc C SMMIMM = sum(I**M) for M, resolution, acentric/centric C NMMNUM = number C SMMTOT = SMMIMM for all resolution bins C NMMTNM = number C THMMTA = theoretical values (acentric) C THMMTC = theoretical values (centric) INTEGER MAXMMT PARAMETER (MAXMMT=4) DOUBLE PRECISION SMMIMM(MAXMMT,KNBINS,2), SMMTOT(MAXMMT), $ SIMEAN,SMMEMM(MAXMMT,KNBINS,2),SMMTOTE(MAXMMT) REAL THMMTA(2:MAXMMT),THMMTC(2:MAXMMT) INTEGER NMMNUM(KNBINS,2), NMMTNM C----- C C Parity groups C MAXPGP number of parity groups C MAXPGC number of parity condition vectors INTEGER MAXPGP,MAXPGC PARAMETER (MAXPGP=8, MAXPGC = MAXPGP + 2) REAL SUMIPG(2,MAXPGP,KNBINS),EOMEAN(2,MAXPGP),EOSUM(2,MAXPGP) INTEGER NRFPGP(2,MAXPGP,KNBINS),NRFSUM(2,MAXPGP),NNNPGP INTEGER HMIN,HMAX,KMIN,KMAX,LMIN,LMAX,IERR C C Groups are defined by a vector, so that the parity is C MOD(scalar_product(hkl.vector), vector(4,)) REAL LSTLSQ EXTERNAL LSTLSQ INTEGER PGPVEC(4,MAXPGC) CHARACTER*12 LABPGP(MAXPGP) C-harvest INTEGER NDATASETS,ISETS(MAXSETS),CSETID(NLPRGI) CHARACTER PNAME_MTZIN(MAXSETS)*64,DNAME_MTZIN(MAXSETS)*64 CHARACTER PNAME_OUT(NCOLS)*64,DNAME_OUT(NCOLS)*64 CHARACTER PNAME_KEYWRD*64,DNAME_KEYWRD*64 CHARACTER Hversion*30 COMMON /CHARV/ PNAME_KEYWRD,DNAME_KEYWRD LOGICAL PRIVATE,USECWD,NOHARVEST LOGICAL GOT_PNAME,GOT_DNAME INTEGER NROWLIMIT,IvalueNotDet REAL ValueNotDet COMMON /IHARV/PRIVATE,USECWD,NOHARVEST, + GOT_PNAME,GOT_DNAME, + NROWLIMIT,IvalueNotDet, + ValueNotDet C-harvest C Last 3 go together for F-centre test DATA PGPVEC/ C h = 2n, k = 2n, l = 2n $ 1,0,0, 2, 0,1,0, 2, 0,0,1, 2, C h+k = 2n, h+l = 2n, k+l = 2n $ 1,1,0, 2, 1,0,1, 2, 0,1,1, 2, C h+k+l = 2n $ 1,1,1, 2, C h+k = 2n & h+l = 2n & k+l = 2n $ 1,1,0, 2, 1,0,1, 2, 0,1,1, 2/ C DATA LABPGP/ $ ' h ', $ ' k ', $ ' l ', $ ' h+k ', $ ' h+l ', $ ' k+l ', $ ' h+k+l ', $ ' h+k,h+l,k+l'/ C DATA NREFWL/0/,WILSC/1.0,1.0/,WILBT/0.0,0.0/ DATA THMMTA/2.0,6.0,24.0/, THMMTC/3.0,15.0,105.0/ DATA THA/9.5,18.1,25.9,33.0,39.3,45.1,50.3,55.1,59.3,63.2/ DATA THC/24.8,34.5,41.6,47.3,52.1,56.1,59.7,62.9,65.7,68.3/ DATA FF/KNBINS*0.0/,SD/KNBINS*0.0/,N/KNBINS*0/,NW/KNBNS2*0/, + FRACS/KNBS20*0.0/,NTOT/2*0/,FRACT/20*0.0/ C C C---- MTZLOK array containing index from program labels C to file labels. On input these should be C set to -1 for compulsory labels and 1 C for optional labels, and exit labels C which are not present in the file have C their lookup entry set to 0. DATA MTZLOK/3*-1,13*1/ DATA MTZERR/0/,MTZPRT/1/,MTZFLG/0/,MTITOT/0/,LABFLG/0/ DATA PI /3.14159265/ C C---- LSPRGI array containing the program label strings C labels input C DATA (LSPRGI(J),J=1,16)/'H','K','L','IMEAN', + 'SIGIMEAN','IW(+)','SIGIW(+)','I(+)','SIGI(+)', + 'IW(-)','SIGIW(-)','I(-)','SIGI(-)','F','SIGF','FreeR_flag'/ C C---- CTPRGI array containing the program label types C DATA (CTPRGI(J),J=1,16)/'H','H','H','J','Q','K','M','K', + 'M','K','M','K','M','F','Q','I'/ C C---- o/p label types C DATA (LSPRGO(J),J=1,19)/'H','K','L','F','SIGF','DANO','SIGDANO', + 'F(+)','SIGF(+)','F(-)','SIGF(-)', + 'IMEAN','SIGIMEAN','I(+)','SIGI(+)','I(-)','SIGI(-)', + 'FreeR_flag','ISYM'/ DATA (CTPRGO(J),J=1,19)/'H','H','H','F','Q', 'D', 'Q', + 'G', 'L', 'G', 'L', + 'J', 'Q', 'K', 'M', 'K', 'M', + 'I', 'Y'/ C IFAIL = 0 LUNIN = 5 LUNOUT = 6 C C **************************************** CALL CCPFYP CALL MTZINI call ccp4h_init() call ccp4h_header('TRUNCATE','trn',2) call ccp4h_pre_beg() CALL CCPDAT(DATEMT) CALL CCPDPN(LUNIN,'DATA','READONLY','F',LDUM,IFAIL) CALL CCPDPN(LUNOUT,'PRINTER','PRINTER','F',LDUM,IFAIL) C CALL CCPRCS(LUNOUT,'TRUNCATE','1994/10/14 14:47:05') CALL CCPRCS(LUNOUT,'TRUNCATE','$Date: 2002/08/20 14:51:52 $') C **************************************** C IBATCH = 1 LNOTRN = .FALSE. RNGSET = .FALSE. LTRAP = .TRUE. C LOGOUT controls whether or not HKLOUT is output LOGOUT = .TRUE. C C WRITE (LUNOUT,FMT=6000) 6000 FORMAT (5X,'TRUNCATE INTENSITIES TO AMPLITUDES',/, + 5X,'==================================',////) call ccp4h_pre_end() call ccp4h_summary_beg() call ccp4h_toc_beg() call ccp4h_toc_ent('Command Input','#command') call ccp4h_toc_ent('Input MTZ File','#inp') call ccp4h_toc_ent('Output File','#out') call ccp4h_summary_end() call ccp4h_toc_ent('Volume, Solvent Content etc','#volume') call ccp4h_summary_beg() call ccp4h_toc_ent('Scale from Wilson Plot','#wilson') call ccp4h_summary_end() call ccp4h_toc_ent('Analysis of Mean Intenity','#mi') call ccp4h_toc_ent('Header Information to Output MTZ File', + '#header') call ccp4h_toc_ent('Distribution of Observed Intensity', + '#obsd') call ccp4h_summary_beg() call ccp4h_toc_ent('Acentric Moments of Intensity', + '#acentric') call ccp4h_toc_ent('Centric Moments of intensity', + '#centric') call ccp4h_summary_end() call ccp4h_toc_ent('Cumulative Intensity Distribution', + '#cumm') call ccp4h_toc_ent('Mean Amplitude vs. Resolution', + '#resolution') call ccp4h_toc_ent('Anisotropic Analysis: FALLOFF','#trun-anis') call ccp4h_summary_beg() call ccp4h_toc_end() call ccp4h_rule() call ccp4h_header('Command Input','command',3) call ccp4h_pre_beg() call ccp4h_link('TITLE ','truncate.html#title') call ccp4h_link('TRUNCATE ','truncate.html#truncate') call ccp4h_link('NRESIDUE ','truncate.html#nresidue') call ccp4h_link('LABOUT ','truncate.html#labout') call ccp4h_link('ANOMALOUS','truncate.html#anomalous') call ccp4h_link('CELL ','truncate.html#cell') call ccp4h_link('CONTENTS ','truncate.html#contents') call ccp4h_link('HEADER ','truncate.html#header') call ccp4h_link('FALLOFF ','truncate.html#falloff') call ccp4h_link('HISTORY ','truncate.html#history') call ccp4h_link('LABIN ','truncate.html#labin') call ccp4h_link('PLOT ','truncate.html#plot') call ccp4h_link('RANGES ','truncate.html#ranges') call ccp4h_link('RESOLUTION','truncate.html#resolution') call ccp4h_link('RSCALE ','truncate.html#rscale') call ccp4h_link('SCALE ','truncate.html#scale') call ccp4h_link('SYMMETRY ','truncate.html#symmetry') call ccp4h_link('VPAT ','truncate.html#vpat') call ccp4h_summary_end() C C---- Read control cards C C ***** CALL KEYIN C ***** MINBIN = 10000 RMAXI = 0.0 SMIN = 1.0 SMAX = 0.0 FMIN = 1000000.0 FMAX = 0.0 NSYSABS = 0 NREFWL = 0 MTZIN = 1 MTZOUT = 1 MTZERR = 0 DO 5 I=1,KNBINS FLOWER(I) = 0.0 FF(I) = 0.0 FFD2(I) = 0.0 SD(I) = 0.0 N(I) = 0 NRFUNQ(I) = 0 NPOSS(I) = 0 SUMIIS(I) = 0.0 SUMISS(I) = 0.0 SUMIRS(I) = 0.0 SUMINS(I) = 0.0 SUMFIS(I) = 0.0 SUMFSS(I) = 0.0 SUMFRS(I) = 0.0 DO 6, J=1,MAXPGP SUMIPG(1,J,I) = 0.0 SUMIPG(2,J,I) = 0.0 NRFPGP(1,J,I) = 0 NRFPGP(2,J,I) = 0 6 CONTINUE NMMNUM(I,1) = 0 NMMNUM(I,2) = 0 DO 7, J = 1,MAXMMT SMMEMM(J,I,1) = 0.0 SMMEMM(J,I,2) = 0.0 SMMIMM(J,I,1) = 0.0 SMMIMM(J,I,2) = 0.0 7 CONTINUE 5 CONTINUE DO 8, J=1,MAXPGP EOSUM(1,J) = 0.0 EOSUM(2,J) = 0.0 8 CONTINUE call ccp4h_pre_end() call ccp4h_rule() call ccp4h_summary_beg() call ccp4h_header('Input MTZ File','inp',3) call ccp4h_summary_end() call ccp4h_pre_beg() C C---- Open input MTZ dataset C C HEADER "string" Controls printout from reading MTZ HEADER C string is one of the following; C NONE sets MTZPRT = 0 (default) C no header o/p C BRIEF sets MTZPRT = 1 C brief header o/p C HIST sets MTZPRT = 2 C brief + mtz history C ALL sets MTZPRT = 3 C full header o/p from mtz reads C C ======================================== CALL LROPEN(MTZIN,'HKLIN',MTZPRT,MTZERR) C ======================================== C IF (MTZERR.EQ.-1) THEN WRITE (LUNOUT,FMT=*) ' LROPEN no such file for HKLIN' CALL CCPERR(1,' LROPEN no such file for HKLIN') END IF C C---- CALL to return information about the MTZ file open for C read on index INDEX. C C ============================================== CALL LRINFO(MTZIN,VERSNX,NUMCOL,NREFLX,RNGMTZ) C ============================================== C C WRITE (LUNOUT,FMT=6002) VERSNX,NUMCOL,NREFLX C 6002 FORMAT (/,' Reading from HKLIN mtz file:',/, C + ' This file written with MTZLIB Version Number : ',A,/, C + ' File HKLIN contains a total of ',I6,' Columns',/, C + ' and a total of ',I10,' Reflections') C C---- CALL to return the title and it's length from the C header of an MTZ file. C C ========================== CALL LRTITL(MTZIN,TITIN,LENTIT) C ========================== C IF (LENTIT.LE.1) TITIN = ' From Truncate on the '//DATEMT C C---- CALL to READ Cell Parameters from the header common block C C ==================== CALL LRCELL(MTZIN,CELMTZ) C ==================== C C---- Call to read resolution limits, if not read from commands IF (.NOT.KYRESO) THEN CALL LRRSOL(MTZIN,SMINWL,SMAXWL) RMINWL = 10000 RMAXWL = 0.5 IF (SMINWL.GT.0.0) RMINWL=1./SQRT(SMINWL) IF (SMAXWL.GT.0.0) RMAXWL=1./SQRT(SMAXWL) C WRITE(LUNOUT,'(A,4F10.4)') C + ' Input File resolution limits rminwl ..sminwl..', C + RMINWL,RMAXWL,SMINWL,SMAXWL ENDIF C C---- Check if user wants to override HKLIN symmetry with C Key_Word input symmetry C IF (KYSYMM) THEN C C---- Set up symmetry (RDSYMM must be called first ) C C ********************************************************* CALL ASUSET(NAMSPG,NSPGRP,PGNAME,NSYM,RSYM,NSYMP,NLAUE,.TRUE.) C ********************************************************* C ELSE C C---- CALL to return symmeytry information (other than symmetry C operations) from the MTZ header C C ==================================================== CALL LRSYMI(MTZIN,NSYMP,LATTYP,NSPGRP,NAMSPG,PGNAME) C ==================================================== C C IF (NSPGRP.GT.0) THEN C C C---- CALL to return symmetry operations from header of MTZ file C C ============================ CALL LRSYMM(MTZIN,NSYM,RSYM) C ============================ C C ************************************************************** CALL ASUSET(NAMSPG,NSPGRP,PGNAME,NSYM,RSYM,NSYMP,NLAUE,.TRUE.) C ************************************************************** C ELSE WRITE (LUNOUT,FMT=6004) 6004 FORMAT (/,' NO SYMMETRY in HKLIN and NO Key_Word input ', + 'for SYMMETRY ... STOPPING !!!') CALL CCPERR(1,'No symmetry') END IF END IF C C---- Now I can work out centricity and epsilon classes C This will replace icent(...) and Centric_zone C C C ***************************** CALL CENTRIC(NSYM,RSYM,0) CALL EPSLN(NSYM,NSYMP,RSYM,0) C ***************************** C C---- Set up reciprocal parameters for THETA calculation C Test if user wants to over ride HKLIN cell values C with Key_Word CELL C IF (CELL(1).LE.0.001) THEN C C DO 41 JDO = 1,6 CELL(JDO) = CELMTZ(JDO) 41 CONTINUE C C END IF C C ***************** CALL SETCLL(CELL) C ***************** C WRITE (LUNOUT,FMT=6006) CELL 6006 FORMAT (/,6X,'Cell Dimensions: ',6F12.2,/) C Check for possible twinning KYTWIN = .FALSE. IF(PGNAME(1:3).EQ.'PG3' .OR. PGNAME(1:3).EQ.'PG4' $.OR. PGNAME(1:3).EQ.'PG6' .OR. PGNAME(1:4).EQ.'PG23') $ KYTWIN = .TRUE. IF(PGNAME(1:3).EQ.'PG2' ) THEN C If cos (betast) = n|c|/2|a| or |a|/(2n|c|) twinning possible C Either |a+nc| = |a| or |a+nc| = |nc| C Test n=1,2. Value of 0.015 a bit arbitrary! COSBST = -cos(PI*CELL(5)/180.0) IF (ABS(COSBST-0.50*CELL(3)/CELL(1)).LT.0.015 $ .OR. ABS(COSBST- CELL(3)/CELL(1)).LT.0.015 $ .OR. ABS(COSBST-0.50*CELL(1)/CELL(3)).LT.0.015 $ .OR. ABS(COSBST-0.25*CELL(1)/CELL(3)).LT.0.015 ) KYTWIN = .TRUE. END IF if(kytwin) CALL CCPERR(2, + ' **** Beware! - Cell dimensions could permit Twinning ****') C FALLOFF: set variables for s/r TRICART IF (KYFALL) THEN aa=cell(1) bb=cell(2) cc=cell(3) alpha=cell(4) beta=cell(5) gamma=cell(6) rad=pi/180. alpha=alpha*rad beta=beta*rad gamma=gamma*rad fvol=aa*bb*cc*sqrt(1.0-(cos(alpha))**2-(cos(beta))**2- 1 (cos(gamma))**2+2.0*cos(alpha)*cos(beta)*cos(gamma)) ast=(bb*cc*sin(alpha))/fvol bst=(aa*cc*sin(beta))/fvol cst=(aa*bb*sin(gamma))/fvol alphast=acos((cos(beta)*cos(gamma)-cos(alpha)) 1 /(sin(beta)*sin(gamma))) betast=acos((cos(alpha)*cos(gamma)-cos(beta)) 1 /(sin(alpha)*sin(gamma))) gammast=acos((cos(alpha)*cos(beta)-cos(gamma)) 1 /(sin(alpha)*sin(beta))) c---- c----calculate the matrix that transforms the cell coordinates h,k,l to c----Cartesian coordinates. c---- call tricart (ast,bst,cst,alphast,betast,gammast,transf) c---- c----Set up resolution bins and initialize arrays c---- sthmax=.25*SMAXWL sthmin=.25*SMINWL sthstap=(sthmax)/REAL(NRANGE) nmax=NRANGE nmin=int(sthmin/sthstap)+1 if(nmax.gt.KNBINS) + CALL CCPERR(1,' Array overflow (nmax.gt.KNBINS) ') do 97 nsinth=1,KNBINS somov(nsinth)=0.0 somsdov(nsinth)=0.0 numov(nsinth)=0 97 continue do 98 j=1,3 do 99 nsinth=1,KNBINS somdir(j,nsinth)=0.0 somsddir(j,nsinth)=0.0 numdir(j,nsinth)=0 99 continue 98 continue ENDIF C C---- CALL to apply the input column label assignments to the C MTZ file which is open for read and setup the Lookup pointer array C Also note that having an unassigned program label is a C fatal error, but this CALL uses LERROR to report C this as a warning, so that all the program label assignments C can be checked before stopping. C C ============================================== CALL LRASSN(MTZIN,LSPRGI,NLPRGI,MTZLOK,CTPRGI) C ============================================== C If amplitude input no output generated; bad idea to TRUNCATE a set of Fs C Also, don't apply Wilson scale. SCALEF = 100 unless given on SCALE card. IF(MTZLOK(14).NE.0) THEN LOGOUT = .FALSE. LAPSCL = .FALSE. IF (SCALEF .EQ. 0.0) SCALEF = 100.0 SCALEFINV = 1.0/(SCALEF*SCALEF) ENDIF C---- LSPRGI array containing the program label strings C labels input C C DATA (LSPRGI(J),J=1,16)/'H','K','L','IMEAN', C + 'SIGIMEAN','IW(+)','SIGIW(+)','I(+)','SIGI(+)', C + 'IW(-)','SIGIW(-)','I(-)','SIGI(-)','F','SIGF','FreeR_flag'/ C C C C C C---- Check that data pairs are either both present or both absent C IERR = 0 DO I = 6,14,2 NCIP = MTZLOK(I) NCSP = MTZLOK(I+1) IF ((NCIP+NCSP).NE.0 .AND. (NCIP*NCSP).EQ.0) THEN WRITE (LUNOUT,FMT=6010) LSPRGI(I),LSPRGI(I+1) 6010 FORMAT (//, $ ' !!!! Columns ',A,' and ',A,' must be either both pres', + 'ent or both absent !!!!',/) IERR = 1 ENDIF END DO IF(IERR.EQ.1) CALL CCPERR(1,'Unmatched data pair') C C---- Check if IW+/- and/or I+/- columns are present in input file. C Default to include anomalous if present on input file C IF (KANOM .EQ. -1) THEN KANOM = MTZLOK(6)*MTZLOK(7)*MTZLOK(10)*MTZLOK(11) + + MTZLOK(8)*MTZLOK(9)*MTZLOK(12)*MTZLOK(13) ENDIF IF (KANOM .LE. 0) THEN c DATA (LSPRGO(J),J=1,19)/'H','K','L','F','SIGF','DANO','SIGDANO', c + 'F(+)','SIGF(+)','F(-)','SIGF(-)', c + 'IMEAN','SIGIMEAN','I(+)','SIGI(+)','I(-)','SIGI(-)', c + 'FreeR_flag','ISYM'/ C---- No anomalous data present; rearrange output columns LSPRGO(6)=LSPRGO(12) CTPRGO(6)=CTPRGO(12) LSPRGO(7)=LSPRGO(13) CTPRGO(7)=CTPRGO(13) NLPRGO = 7 JLOOP = 1 ELSE C---- Anomalous data present NLPRGO = 18 JLOOP = 9 END IF C C#### Copy FreeR_flag column if present. IF (MTZLOK(16).EQ.0) THEN IF (KANOM.GT.0) THEN LSPRGO(18)=LSPRGO(19) CTPRGO(18)=CTPRGO(19) ENDIF ELSE NLPRGO = NLPRGO + 1 IF (KANOM.LE.0) THEN LSPRGO(8) = LSPRGO(18) CTPRGO(8) = CTPRGO(18) ENDIF ENDIF C call ccp4h_rule() call ccp4h_pre_end() call ccp4h_summary_beg() call ccp4h_header('Output MTZ File','out',3) IF (.NOT.LOGOUT) WRITE(6,'(A)') + ' No output MTZ file written. ' call ccp4h_summary_end() call ccp4h_pre_beg() C C---- Open output MTZ file C C *********************** IF(LOGOUT) CALL LWOPEN(MTZOUT,'HKLOUT') C *********************** C C---- CALL to update the title of an MTZ file in the header C common block C C C TITIN from HKLIN file C TITOUT from Key_Word file C C---- AND C C---- CALL to write new history lines to history header. C History headers have a maximum of NHISTL lines, we write the new C NLINES of history to the header and then fill up any free lines C with the older lines. C C C---- If TITLE Key_Word used then put in New Title C IF (.NOT.KYTITL) TITOUT = TITIN C C ============================ IF(LOGOUT) CALL LWTITL(MTZOUT,TITOUT,MTITOT) C ============================ C C---- IF HISTORY Key_Word used then add new history to old C IF (.NOT.KYHIST) THEN C C---- No Key_Word HISTORY therefore just add one new history line C IF(LOGOUT) CALL LWHSTL (MTZOUT, ' ') ELSE IF(LOGOUT) CALL LWHIST(MTZOUT,HISOUT,NHSOUT) END IF C C---- CALL to write Cell Parameters into the header common block C Only write explicitly if changed on input C =================== IF (KYCELL.AND.LOGOUT) CALL LWCELL(MTZOUT,CELL) C =================== C C---- CALL to update the symmetry operations and information C in the MTZ header (only if symmetry read from input) C LATTYP = NAMSPG(1:1) C C ========================================================= IF (KYSYMM.AND.LOGOUT) + CALL LWSYMM(MTZOUT,NSYM,NSYMP,RSYM,LATTYP,NSPGRP,NAMSPG,PGNAME) C ========================================================= C c-harvest C C---- If PNAME and DNAME cards were given, get new dataset names C IF (GOT_PNAME .AND. GOT_DNAME) THEN C C---- ok can go ahead NCOLS set but only NLPRGO used C DO 115 Jcol = 1,NCOLS PNAME_OUT(Jcol) = PNAME_KEYWRD DNAME_OUT(Jcol) = DNAME_KEYWRD 115 CONTINUE ELSE C C---- must have PNAME and DNAME from MTZ header as C Read in list of datasets in input file C Assign everything to 1st dataset; shouldn't have more than one! C CALL LRID(MTZIN,PNAME_MTZIN,DNAME_MTZIN,ISETS,NDATASETS) IF (NDATASETS.GT.0) THEN CALL LRCLID(MTZIN,CSETID,NUMCOL) C C---- Assign protein name/dataset name for output columns C Cycle through program labels that are assigned C DO 3115 Jcol = 1,NCOLS PNAME_OUT(Jcol) = PNAME_MTZIN(1) DNAME_OUT(Jcol) = DNAME_MTZIN(1) 3115 CONTINUE C---- Single pname/dname used for harvesting PNAME_KEYWRD = PNAME_MTZIN(1) DNAME_KEYWRD = DNAME_MTZIN(1) ENDIF END IF C C---- Now want column labels and column numbers for C output file IF(LOGOUT) + CALL LKYOUT(MTZOUT,LSPRGO,NLPRGO,LBOTKN,LBOLIN,LBOBEG,LBOEND) IAPPND = 0 IF(LOGOUT)CALL LWASSN(MTZOUT,LSPRGO,NLPRGO,CTPRGO,IAPPND) IF (PNAME_KEYWRD.NE.' ') THEN C C---- Add in new dataset. Library will remove duplicates. C---- Store the project name and dataset name in the mtz header: C IF(LOGOUT)CALL LWID(MTZOUT,PNAME_KEYWRD,DNAME_KEYWRD) IF(LOGOUT)CALL LWIDAS(MTZOUT,NLPRGO,PNAME_OUT,DNAME_OUT,IAPPND) ENDIF C C the ccp4 mtz nan for missing data is not very good C harvest will need to check for not-determined values C and for calculated values that are not 'finite' CALL CCP4_VERSION(Hversion) Hversion = 'CCP4_' // Hversion(1:lenstr(Hversion)) C If NOHARVEST keyword specified, switch off harvesting by setting C ProjectName and DataSetName to ' ' C Also, no harvesting if we're just analysing existing structure factors. IF (NOHARVEST .OR. (.NOT.LOGOUT)) THEN PNAME_KEYWRD = ' ' DNAME_KEYWRD = ' ' ENDIF C C---- HarvestInit C C SoftwareClass = description of program function C Package = CCP4 C ProgramName = mlphare C ProgramVersion = 'CCP4_3.1 of 1996/11/01' C ProjectName = Name of protein/complex under study C DataSetName = Name of key data set for which cell etc apply to C useCWD if true write to current working directory C Private if true chmod DepositFiles directory "drwx------" C IvalueNotDet RETURNED value to which int variables C are initialised C and if unchanged no data_item written C ValueNotDet RETURNED value to which real variables C are initialised C and if unchanged no data_item written C NRowLimit either 80 or 132 is max width of a row on o/p C C writes C C _entry.id ProjectName C _entry.data_set_id DataSetName C _audit.creation_DayTime as fortran fdate string C _software.name ProgramName C _software.version ProgramVersion C _software.classification SoftwareClass C C If useCWD false then C if directory $HOME/DepositFiles doesnt exist it will be created C with protection drwxr-xr-x C or if Private set drwx----- C if directory $HOME/DepositFiles/ProjectName doesnt exist it will C be created C harvest file then written to C $HOME/DepositFiles/ProjectName/DataSetName.ProgramName C overwritting any file that exists with this name C else if useCWD true C the harvest file written to C ./DataSetName.ProgramName C C Package is used internally for CCP4 style of NAN C will be used for other package specific rules C C IvalueNotDet and ValueNotDet used to decide if a data_item C is to be written C C Internal check done on if a calculated value is finite C if not no data_item written C C---- if no PROJECTNAME and/or DATASETNAME given then look in MTZ header C CALL Hinitialise( + 'CCP4', + 'truncate', + Hversion, + PNAME_KEYWRD,DNAME_KEYWRD, + useCWD, + Private, + IvalueNotDet, + ValueNotDet, + NRowLimit) C C _software.classification SoftwareClass C _item_enumeration.value 'data collection' C 'data reduction' C 'phasing' C 'model building' C 'refinement' C 'validation' C 'other' C _software.contact_author C _software.contact_author_email C _software.description C CALL Hsoftware('data reduction','K.S. Wilson and S. French', + 'ccp4@dl.ac.uk', +'calculates a best estimate of scaled F from I & sd(I)') C C---- if a standard symmetry setting -only need IntTabNum C else send SpaceGrpNam, number_equivalent_positions, C and symmetry operations, and in this case IntTabNum C must be greater than 230 C C if non-standard number is given AND it is a known ccp4 C number then you dont need any more info, current C known non-standard ones are C 1003, 3004, 1004, C 2005, 1005, 1018, C 1020, 1021, 1022, C 1023, 1059, 1094, C 1197 C C _Symmetry.Int_Tables_number C _Symmetry.space_group_name_H-M C C loop_ C _Symmetry_equiv.id C _Symmetry_equiv.pos_as_xyz C CALL HSymmetry(NSPGRP,NAMSPG,NSYM,RSYM) C C---- harvest C C _cell.length_a C _cell.length_b C _cell.length_c C _cell.angle_alpha C _cell.angle_beta C _cell.angle_gamma C CALL HCell(cell) C C----harvest C C----- Set up resolution limits if required C RMINWL, RMAXWL, SMINWL SMAXWL should have been defined either by C RESOLUTION command, or from input file C Set limits for scaling if not set explicitly IF (.NOT.SCBLOG) THEN C -- no RSCALE defined, set scaling limits == overall limits RMAXSC = RMAXWL RMINSC = RMINWL IF(RMAXSC.LE.3.5) RMINSC = 4.0 SMAXSC = SMAXWL SMINSC = SMINWL IF(RMAXSC.LE.3.5) SMINSC = 0.0625 ELSE C -- force scaling limits to lie within overall resolution limits SMINSC = MAX(SMINSC,SMINWL) SMAXSC = MIN(SMAXSC,SMAXWL) RMINSC = 1./SQRT(SMINSC) RMAXSC = 1./SQRT(SMAXSC) ENDIF C Check for RMINSC gt than RMAXSC.. IF(RMINSC.LT.RMAXSC) THEN WRITE(6,'(A,A,2F8.2,/,A)') + ' *** The automatic scaling resolution choice is confused -', + ' RMINSC RMAXSC',RMINSC,RMAXSC, + ' Minimum greater than Maximum!! - Define RSCALE yourself **' CALL CCPERR(1,' Fatal error') END IF C C..Set bin width or number of bins. IF (RANGE .LT. 0.0) THEN RANGE = (SMAXWL-SMINWL)/NRANGE ELSE NRANGE = ((SMAXWL-SMINWL)/RANGE)+1 C..If width of bin was set explicitely on input this will often result C in a tail-end of few reflections in the top bin. To avoid this the width C of bin is adjusted so the resolution range is covered evenly. RANGE = (SMAXWL-SMINWL)/NRANGE ENDIF C WRITE(LUNOUT,FMT='(/,A,F7.4,/,A,I4)') + ' Width of bin : ',RANGE,' Number of bins : ',NRANGE C IF (NRANGE.GT.KNBINS) THEN WRITE (LUNOUT,FMT=6762) KNBINS 6762 FORMAT (//,1X,'***** FATAL ERROR *****',/,1X, + 'Reflection outside the ',I4,'th resolution bin...',/,1X, + 'This can result in reflections being incorrectly rejected',/, + 1X,'RERUN with a larger value specified on keyword "RANGES" ', + 'or use default.') CALL CCPERR(1,'REFLECTION(S) OUT OF RESOLUTION RANGE') ENDIF C C..is sampling too coarse? IF (RANGE.GT.WDTHMX) THEN WRITE(LUNOUT,FMT='(A,F7.3,A,/,A)') + ' Width of bin too large (gt ',WDTHMX,').', + ' Rerun the program with more bins or try the default value.' CALL CCPERR(3,'SAMPLING MAY BE TOO COARSE') ENDIF C IF (RANGE.GT.0.02) + WRITE(LUNOUT,'('' WARNING: WIDTH OF BIN LARGE (GT 0.02)'')') C C----------------------- C Check the number of reflections per bin are sensible. ( taken from UNIQUE) C ****************************************** CALL HKLLIM(NLAUE,RMINWL,RMAXWL,CELL,HMIN,HMAX,KMIN,KMAX,LMIN, + LMAX) C ****************************************** C WRITE (6,FMT=6002) HMIN,HMAX,KMIN,KMAX,LMIN,LMAX 6002 FORMAT (/,1X,'Limits on H,K,L..',3 (I4,' to',I4,3X),/) C C 30 CONTINUE C C---- need to initialise these each time NIN = 0 NREJS = 0 NOUT = 0 DO 77 I=1,KNBINS NRFUNQ(I) = 0 NPOSS(I) = 0 77 CONTINUE DO 70 IH = HMIN,HMAX DO 60 IK = KMIN,KMAX DO 50 IL = LMIN,LMAX NIN = NIN + 1 C C---- Test for limiting conditions for indices C INHKL(1) = IH INHKL(2) = IK INHKL(3) = IL C C IF (IH.EQ.0 .AND. IK.EQ.0 .AND. IL.EQ.0) GO TO 50 C C INASU returns +1 if reflection is in bounds C ********************* IF (INASU(INHKL, NLAUE) .LE. 0) THEN GO TO 50 ENDIF C ********************* C C C---- EPS = multiplicity of this reflection. C ISYSAB = 1 - systematic absence C ISYSAB = 0 C C *********************** CALL EPSLON(INHKL,EPS,ISYSAB) C *********************** C C---- Test for systematic absences IF (ISYSAB .EQ. 1) GO TO 50 C C---- Calculate resolution C C ********************** SINTHL = LSTLSQ(MTZIN,IH,IK,IL)*4.0 C ********************** C IF ((SINTHL.LT.SMINWL) .OR. (SINTHL.GT.SMAXWL)) THEN NREJS = NREJS + 1 GO TO 40 END IF C NT = (SINTHL-SMINWL)/RANGE C..this should take care of the reflection(s) JUST on the upper edge IF(NT.LT.NRANGE) NT = NT + 1 C NPOSS(NT) = NPOSS(NT) + 1 NRFUNQ(NT) = NRFUNQ(NT) + NINT(NSYM/EPS) C C NOUT = NOUT + 1 40 CONTINUE C 50 CONTINUE 60 CONTINUE 70 CONTINUE C C C..Check that there are no less than MINREF reflections in any one bin DO 75 NT = 1,NRANGE MINBIN = MIN(NPOSS(NT),MINBIN) 75 CONTINUE C IF (MINBIN.LT.MINREF) THEN IF (.NOT.RNGSET) THEN C..RANGE was not set explicitely: C..Increase bin width NRANGE=NRANGE-1 RANGE = (SMAXWL-SMINWL)/NRANGE MINBIN = 1000 C WRITE(LUNOUT,FMT='(/,A,I6,A,/,A,F7.4,/,A,I4)') + ' Too few reflections per bin.',MINBIN, + ' Increasing bin width:',' New bin width : ',RANGE, + ' Number of bins: ',NRANGE C C..is sampling too coarse? IF (RANGE.GT.WDTHMX) THEN WRITE(LUNOUT,FMT='(A,F7.3,A)') + ' Width of bin too large (gt ',WDTHMX,').' CALL CCPERR(1,'WIDTH OF BIN TOO LARGE') ENDIF C IF (RANGE.GT.0.02) + WRITE(LUNOUT,'('' WARNING: WIDTH OF BIN LARGE (GT 0.02)'')') C C GO TO 30 C..RANGE was set explicitly. Don't do anything but print a warning ELSE WRITE(LUNOUT,FMT='(/,A,I5,A,/,A,/,A)') + ' WARNING: Less than ',MINREF,' reflections in some bin(s).', + ' This may result in bad statistics.', + ' Try rerunning the program without the RANGES card.' ENDIF ENDIF C C C---- Call WILPRE to estimate scattering C C ***************** CALL WILPRE(NSYM,CELL) C ***************** C C---- 1st pass through MTZ file. C Return here to read all MTZ reflections C C---- CALL to read a reflection record from an MTZ file which C has been opened for read. This returns the record in program C order. C C ============================ 110 CALL LRREFF(MTZIN,S,RECIN,MTZEOF) C ============================ C IF (MTZEOF) GO TO 130 C C ==================== CALL LRREFM(MTZIN,LOGMSS) C ==================== C C---- test consistancy of MNFs C IF (LTRAP) THEN IF (LOGMSS(4) .NEQV. LOGMSS(5)) LTRAP = .FALSE. IF (LOGMSS(6) .NEQV. LOGMSS(7)) LTRAP = .FALSE. IF (LOGMSS(8) .NEQV. LOGMSS(9)) LTRAP = .FALSE. IF (LOGMSS(10) .NEQV. LOGMSS(11)) LTRAP = .FALSE. IF (LOGMSS(12) .NEQV. LOGMSS(13)) LTRAP = .FALSE. IF (LOGMSS(14) .NEQV. LOGMSS(15)) LTRAP = .FALSE. C IF (LOGMSS(16) .NEQV. LOGMSS(17)) LTRAP = .FALSE. C IF (LOGMSS(18) .NEQV. LOGMSS(19)) LTRAP = .FALSE. IF (.NOT. LTRAP) WRITE(LUNOUT,*) . ' *** MNFs are inconsistent; could cause unexpected output ***' ENDIF C C---- resolution limits C IF (S.LT.SMINWL .OR. S.GT.SMAXWL) GO TO 110 C C---- include only IMEAN in this first pass through data file C Read file to get WILSON distribution of Intensities C DO 120 JDO = 1,3 INHKL(JDO) = NINT(RECIN(JDO)) 120 CONTINUE C C need EITHER IMEAN and SIGIMEAN or F and SIGF IF ( ((MTZLOK(4).NE.0 .AND. MTZLOK(5).NE.0) + .AND.(LOGMSS(4).OR.LOGMSS(5))) + .OR.((MTZLOK(14).NE.0 .AND. MTZLOK(15).NE.0) + .AND.(LOGMSS(14).OR.LOGMSS(15)))) + GOTO 110 RMAXI = MAX(RMAXI,RECIN(4)) IF (MTZLOK(14).NE.0 .AND. MTZLOK(15).NE.0) THEN F = RECIN(14)*RECIN(14) SS = 2.0*RECIN(15)*RECIN(14) + RECIN(15)*RECIN(15) C Need to reduce scale F = F * SCALEFINV SS = SS * SCALEFINV END IF IF (MTZLOK(4).NE.0 .AND. MTZLOK(5).NE.0) THEN F = RECIN(4) SS = RECIN(5) END IF IF (SS.EQ.0.0) GOTO 110 NT = (S-SMINWL)/RANGE C..this should take care of the reflection(s) JUST on the upper edge IF(NT.LT.NRANGE) NT = NT + 1 C C SMIN = MIN(SMIN,S) SMAX = MAX(SMAX,S) C C---- Compute multiplicity of term for hemisphere of data C C *************************** CALL EPSLON(INHKL,WEIGHT,ISYSAB) C *************************** C FF(NT) = FF(NT) + F/WEIGHT SD(NT) = SD(NT) + SS/WEIGHT N(NT) = N(NT) + 1 AMULT = NSYM/WEIGHT C C---- Look up scattering factor tables generated by ATREC C IBIN = INT(S/DELDSS) + 1 DELIS = S - REAL(IBIN)*DELDSS C C---- Wilson's expected intensity C IF ( IBIN.LE.1) THEN FFscat = FFS(1) ELSE BF=0.5*(FFS(IBIN+1) - FFS(IBIN-1)) AF=BF + FFS(IBIN-1) - FFS(IBIN) FFscat = (AF*DELIS + BF)*DELIS + FFS(IBIN) ENDIF C C Accumulate sums for Wilson plot C SN(NT) = SN(NT) + AMULT SW(NT) = SW(NT) + AMULT*FFscat SR(NT) = SR(NT) + AMULT*S SI(NT) = SI(NT) + AMULT*F C C Accumulate sums for parity groups DO 125, I = 1,MAXPGP C Check parity with respect to this condition IF (I .EQ. MAXPGP) THEN C Special for last one, must pass all 3 tests IS = 2 DO 126, J = I,I+2 K = IDOT(PGPVEC(1,J),INHKL) IF (MOD(K,PGPVEC(4,J)) .NE. 0) GO TO 127 126 CONTINUE IS = 1 ELSE K = IDOT(PGPVEC(1,I),INHKL) IF (MOD(K,PGPVEC(4,I)) .EQ. 0) THEN IS = 1 ELSE IS = 2 ENDIF ENDIF 127 CONTINUE SUMIPG(IS,I,NT) = SUMIPG(IS,I,NT) + F/SS NRFPGP(IS,I,NT) = NRFPGP(IS,I,NT) + 1 125 CONTINUE C C GO TO 110 C C---- End of file reached on 1st pass C 130 CONTINUE C C..Check that there are no less than MINREF reflections in any one bin NREFWL = 0 DO 387 NT = 1,NRANGE MINBIN = MIN(N(NT),MINBIN) NREFWL = NREFWL + N(NT) 387 CONTINUE C C---- Call WILPLT to do Wilson plot C C ************************************************ CALL WILPLT(NSYM,NSYMP,NREFWL,NRANGE,WILSC,WILBT, + IWITCH) C ************************************************ C---- Means of Intensity in radial ranges . C DO 140 NT = 1,NRANGE IF (N(NT).NE.0) THEN FN = REAL(N(NT)) IF(N(NT).GT.0) FF(NT) = FF(NT)/FN IF(N(NT).GT.0) SD(NT) = SD(NT)/FN END IF FLOWER(NT) = REAL(NT-1)*RANGE + SMINWL 140 CONTINUE C C..midpoint of the bins DO 347 NT=1,NRANGE RNGCEN(NT)=FLOWER(NT)+RANGE/2. 347 CONTINUE C C..calculate second derivatives for the spline CALL SCNDER(RNGCEN,FF,NRANGE,FFD2) call ccp4h_rule() call ccp4h_header('Analysis of Mean Intensity','mi',3) call ccp4h_pre_beg() C C---- Output means to line printer . C WRITE (LUNOUT,FMT=6014) 6014 FORMAT (/,' Range Min. S Max. S Dmax(A) Mn(I)/w Mn(SD)', + ' Nref Nposs',/) C C DO 170 NT = 1,NRANGE IF (N(NT).NE.0) THEN UPPER = FLOWER(NT) + RANGE DMAX = 1.0/SQRT(UPPER) WRITE (LUNOUT,FMT=6016) NT,FLOWER(NT),UPPER,DMAX,FF(NT), + SD(NT),N(NT),NPOSS(NT) 6016 FORMAT (1X,I5,2F10.5,F10.2,F10.1,F7.1,2I9) END IF 170 CONTINUE C C---- Output parity group analysis C WRITE (LUNOUT, FMT=6017) (I,I=1,MAXPGP),(LABPGP(I),I=1,MAXPGP) 6017 FORMAT(///, $' Analysis of mean intensity by parity for reflection classes',//, $' For each class, Mn(I/sig(I)) is given for even and odd parity', $' with respect to the condition,',/, $ 'eg group 1: h even & odd; group 7 h+k+l even & odd;', $ ' group 8 h+k=2n & h+l=2n & k+l=2n or not',//, $ ' Range Min_S Dmax Nref',I6,7I12,/, $ 1X,29X,9A,/) C C DO 174, K = 1, MAXPGP NRFSUM(1,K) = 0 NRFSUM(2,K) = 0 174 CONTINUE C DO 171 NT = 1,NRANGE IF (N(NT).NE.0) THEN UPPER = FLOWER(NT) + RANGE DMAX = 1.0/SQRT(UPPER) DO 172, K = 1, MAXPGP DO 175, J = 1,2 IF (NRFPGP(J,K,NT) .GT. 0) THEN EOMEAN(J,K) = SUMIPG(J,K,NT)/NRFPGP(J,K,NT) EOSUM(J,K) = EOSUM(J,K) + SUMIPG(J,K,NT) NRFSUM(J,K) = NRFSUM(J,K) + NRFPGP(J,K,NT) ELSE EOMEAN(J,K) = 0.0 ENDIF 175 CONTINUE 172 CONTINUE WRITE (LUNOUT,FMT=6118) NT,FLOWER(NT),DMAX, + N(NT),((EOMEAN(I,K),I=1,2), K=1,MAXPGP) 6118 FORMAT (1X,I5,F10.5,F7.2,I8,8(2F5.1,2X)) END IF 171 CONTINUE NNNPGP = 0 DO 173, K = 1, MAXPGP NNNPGP = NRFSUM(1,K) + NRFSUM(2,K) EOMEAN(1,K) = 0.0 EOMEAN(2,K) = 0.0 IF (NRFSUM(1,K) .GT. 0) EOMEAN(1,K) = EOSUM(1,K)/NRFSUM(1,K) IF (NRFSUM(2,K) .GT. 0) EOMEAN(2,K) = EOSUM(2,K)/NRFSUM(2,K) 173 CONTINUE WRITE (LUNOUT,FMT=6119) + NNNPGP,((EOMEAN(I,K),I=1,2), K=1,MAXPGP) 6119 FORMAT (/,1X,' Totals: ',I8,8(2F5.1,2X)) C C---- Rewind input file C --- Subroutine to rewind an MTZ file for re-reading of reflections. C The file must already have been opened for read with LROPEN. C The file is postioned at the start of the reflection records. C C C ============ CALL LRREWD(MTZIN) C ============ C NREF = 0 NREJ = 0 IOUT = 0 C C---- Get scale factor for F's from maximum I default = 150. C Amplitudes obtained from the intensities are scaled up by 150.0 C continue using the full range to 32767. C IF (LAPSCL) THEN SCALEF = SQRT(WILSC(1)) ELSE IF (SCALEF .EQ. 0.0) THEN SCALEF = 1.0 ELSE WRITE (LUNOUT, '(/A)') ' Explicit scale given' ENDIF ENDIF SCALEFINV = 1.0/(SCALEF*SCALEF) WRITE (LUNOUT,FMT=6018) SCALEF 6018 FORMAT (/,' Amplitudes will be scaled by',F10.3,' from sqrt(I)',/) C C---- 2nd pass through MTZ file. C Apply truncation and/or sqrt correction. C C ========================== 180 CALL LRREFF(MTZIN,S,RECIN,MTZEOF) C ========================== C IF (MTZEOF) GO TO 230 C C ==================== CALL LRREFM(MTZIN,LOGMSS) C ==================== C C---- resolution limits IF (S.LT.SMINWL .OR. S.GT.SMAXWL) GO TO 180 C C---- Store new reflection, first clear & set arrays/flags. C NREF = NREF + 1 C C IF(LOGOUT)CALL EQUAL_MAGIC(MTZIN,RECOUT,NLPRGO) DO 190 JDO = 1,10 FA(JDO) = 0.0 190 CONTINUE DO 200 JDO = 1,3 INHKL(JDO) = NINT(RECIN(JDO)) IF(LOGOUT)RECOUT(JDO) = RECIN(JDO) 200 CONTINUE C C NT = (S-SMINWL)/RANGE C..this should take care of the reflection(s) JUST on the upper edge IF(NT.LT.NRANGE) NT = NT + 1 IF (NT.GT.NRANGE) NT = NRANGE C C..get the intepolated value CALL SPLINT(RNGCEN,FF,FFD2,NRANGE,S,SIGMA) SIGMA = MAX(SIGMA,0.1) C C---- Test for centric reflection C ICEN = 1 C C ************ CALL CENTR(INHKL,IC) C ************ C IF (IC.EQ.1) ICEN = 2 C C ************************ CALL EPSLON(INHKL,WEIGHT,ISYSAB) C ************************ C SQWT = SQRT(WEIGHT) C C---- Get i and sd C Anomalous data present - JLOOP = 9, 5 pairs of I/sigI C No anomalous data present - JLOOP = 1, 1 pair of I/sigI C DATA (LSPRGI(J),J=1,13)/'H','K','L','IMEAN', C + 'SIGIMEAN','IW(+)','SIGIW(+)','I(+)','SIGI(+)', C + 'IW(-)','SIGIW(-)','I(-)','SIGI(-)','F','SIGF', C + 'F(+)','SIGF(+)','F(-)','SIGF(-)'/ C C Use as a check NFRAC = -1 FFA(1) = 0.0 SSA(1) = 0.0 FFA(2) = 0.0 SSA(2) = 0.0 FFA(3) = 0.0 SSA(3) = 0.0 FFA(4) = 0.0 SSA(4) = 0.0 FFA(5) = 0.0 SSA(5) = 0.0 FA(1) = -1.0 FA(2) = 0.0 FA(3) = 0.0 FA(4) = 0.0 FA(5) = 0.0 FA(6) = 0.0 FA(7) = 0.0 FA(8) = 0.0 FA(9) = 0.0 FA(10) = 0.0 DO 220 JDO220 = 1,JLOOP,2 C Are you inputting F SIGF etc? IF ( JDO220 .EQ.1) THEN IF (MTZLOK(14) .GT.0 .AND..NOT.(LOGMSS(14).OR.LOGMSS(15))) THEN FA(1) = RECIN(14) FA(2) = RECIN(15) FFA(1) = RECIN(14)*RECIN(14)*SCALEFINV SSA(1) = 2.0*RECIN(15)*RECIN(14) + RECIN(15)*RECIN(15) SSA(1) = SCALEFINV*SSA(1) C C---- Totals for Wilson's distribution , Imean only XIA(1) = FFA(1)/WEIGHT C NFRAC = ((XIA(1)*10.0)/SIGMA) NFRAC = NFRAC + 1 IF (NFRAC.LT.1) NFRAC = 1 IF (NFRAC.LE.10) THEN FRACS(NT,NFRAC,ICEN) = FRACS(NT,NFRAC,ICEN) + 1 FRACT(NFRAC,ICEN) = FRACT(NFRAC,ICEN) + 1 END IF NW(NT,ICEN) = NW(NT,ICEN) + 1 NTOT(ICEN) = NTOT(ICEN) + 1 END IF C END IF C inputting Is; if no Is input go straight to 220 IF (MTZLOK(JDO220+3) .EQ.0 .OR. + LOGMSS(JDO220+3) .OR. LOGMSS(JDO220+4)) GOTO 220 FFA(JDO220) = RECIN(JDO220+3) SSA(JDO220) = RECIN(JDO220+4) C IF (SSA(JDO220).LE.0.00000001) GO TO 220 C XIA(JDO220) = FFA(JDO220)/WEIGHT SDIA(JDO220) = SSA(JDO220)/WEIGHT C C---- Totals for Wilson's distribution , Imean only C IF (JDO220.LE.1 .AND. NFRAC.EQ.-1) THEN NFRAC = ((XIA(JDO220)*10.0)/SIGMA) NFRAC = NFRAC + 1 IF (NFRAC.LT.1) NFRAC = 1 IF (NFRAC.LE.10) THEN FRACS(NT,NFRAC,ICEN) = FRACS(NT,NFRAC,ICEN) + 1 FRACT(NFRAC,ICEN) = FRACT(NFRAC,ICEN) + 1 END IF NW(NT,ICEN) = NW(NT,ICEN) + 1 NTOT(ICEN) = NTOT(ICEN) + 1 END IF C C---- Test for TRUNCATION or NOTRUNCATION. C IF (LNOTRN) THEN C C---- NO TRUNCATION option C IF (FFA(JDO220).LT.0.0) FFA(JDO220) = 0.0 C No need to square root if F input FA(JDO220) = SQRT(FFA(JDO220))*SCALEF FA(JDO220+1) = SQRT(SSA(JDO220)+FFA(JDO220))*SCALEF - + FA(JDO220) ELSE IREJ = -1 C C ************************************** IF (ICEN.EQ.1) CALL NEGIAS(XIA(JDO220),SDIA(JDO220),SIGMA,XJ, + SDJ,XF,SDF,IREJ) IF (ICEN.EQ.2) CALL NEGICS(XIA(JDO220),SDIA(JDO220),SIGMA,XJ, + SDJ,XF,SDF,IREJ) C ************************************** C IF (IREJ.EQ.0) THEN FA(JDO220) = XF*SCALEF*SQWT FA(JDO220+1) = SDF*SCALEF*SQWT ELSE C REJECT very negative reflections NREJ = NREJ + 1 FFA(JDO220) = 0.0 GO TO 180 END IF END IF C C 220 CONTINUE C C---- Sums for sd statistics (from IMEAN and SIGIMEAN) IF (SSA(1) .GT. 0.000001) THEN SUMIIS(NT) = SUMIIS(NT) + FFA(1) SUMISS(NT) = SUMISS(NT) + SSA(1) SUMIRS(NT) = SUMIRS(NT) + FFA(1)/SSA(1) SUMINS(NT) = SUMINS(NT) + 1 SUMFIS(NT) = SUMFIS(NT) + FA(1) SUMFSS(NT) = SUMFSS(NT) + FA(2) SUMFRS(NT) = SUMFRS(NT) + FA(1)/FA(2) C---- Sums for moments of I F = FFA(1)/WEIGHT IF(F.GT.0.0)SMMEMM(1,NT,ICEN) = SMMEMM(1,NT,ICEN) + sqrt(F) IF(F.GT.0.0)SMMEMM(3,NT,ICEN) = SMMEMM(3,NT,ICEN) + F*sqrt(F) SMMIMM(1,NT,ICEN) = SMMIMM(1,NT,ICEN) + F SMMIMM(2,NT,ICEN) = SMMIMM(2,NT,ICEN) + F**2 SMMIMM(3,NT,ICEN) = SMMIMM(3,NT,ICEN) + F**3 SMMIMM(4,NT,ICEN) = SMMIMM(4,NT,ICEN) + F**4 NMMNUM(NT,ICEN) = NMMNUM(NT,ICEN) + 1 ENDIF C FALLOFF: IF (KYFALL) THEN sinth=.25*S nsinth=int(sinth/sthstap)+1 IF (NSINTH.LT.NMIN .OR. NSINTH.GT.NMAX) GOTO 211 c---- c----convert h,k,l into Cartesian coordinates. c---- do k=1,3 hc(k)=0.0 do j=1,3 hc(k)=hc(k)+transf(k,j)*REAL(INHKL(j)) enddo enddo IF (FA(1).GE.0.0) THEN do j=1,3 cosang = abs(hc(j))/sqrt(S) C--- cosang can stray just past 1.0 if (cosang.gt.1.0) cosang = 1.0 ang(j)=acos(cosang) if(ang(j).lt.cone*rad) then somdir(j,nsinth)=somdir(j,nsinth)+FA(1) somsddir(j,nsinth)=somsddir(j,nsinth)+FA(1)/FA(2) numdir(j,nsinth)=numdir(j,nsinth)+1 endif enddo somov(nsinth)=somov(nsinth)+FA(1) somsdov(nsinth)=somsdov(nsinth)+FA(1)/FA(2) numov(nsinth)=numov(nsinth)+1 endif 211 continue ENDIF C----Calculate RECOUT columns for HKLOUT if required IF(LOGOUT) THEN C c DATA (LSPRGO(J),J=1,19)/'H','K','L','F','SIGF','DANO','SIGDANO', c + 'F+','SIGF+','F-','SIGF-', c + 'IMEAN','SIGIMEAN','I+','SIGI+','I-','SIGI-','FreeR_flag', c + 'ISYM'/ C C----No anomalous data IF(KANOM.EQ.0) THEN C Only 7 columns output; 6 and 7 are IMEAN and SIGIMEAN IF (.NOT.(LOGMSS(4) .OR. LOGMSS(5))) THEN IF (RECIN(5).NE.0.0) THEN RECOUT(4) = FA(1) RECOUT(5) = FA(2) RECOUT(6) = 0.01*RECIN(4)*SCALEF*SCALEF RECOUT(7) = 0.01*RECIN(5)*SCALEF*SCALEF ENDIF ENDIF RECOUT(8) = RECIN(16) C----Anomalous data present ELSE C Imean and sigImean IF (.NOT.(LOGMSS(4) .OR. LOGMSS(5))) THEN IF (RECIN(5).NE.0.0) THEN RECOUT(12) = 0.01*RECIN(4)*SCALEF*SCALEF RECOUT(13) = 0.01*RECIN(5)*SCALEF*SCALEF ENDIF ENDIF C C F+ and sigF+, I+ and sigI+ IF (.NOT.(LOGMSS(8) .OR. LOGMSS(9))) THEN IF (RECIN(9).NE.0.0) THEN RECOUT(8) = FA(5) RECOUT(9) = FA(6) RECOUT(14) = 0.01*RECIN(8)*SCALEF*SCALEF RECOUT(15) = 0.01*RECIN(9)*SCALEF*SCALEF END IF END IF C C F- and sigF-, I- and sigI- IF (.NOT.(LOGMSS(12).OR.LOGMSS(13))) THEN IF (RECIN(13).NE.0.0) THEN RECOUT(10) = FA(9) RECOUT(11) = FA(10) RECOUT(16) = 0.01*RECIN(12)*SCALEF*SCALEF RECOUT(17) = 0.01*RECIN(13)*SCALEF*SCALEF END IF END IF C C Now want to output Fmean and D .... C Try to derive from I+ and I- (don't bother with IW+ and IW- C which should soon be obsolete). C We assume that at least one of I+ and I- is present. If both C are absent, then presumably so is IMEAN and everything will be MNF. ISYMB = 0 C Case 1: I+ and I- present. IF (FA(6).NE.0.0 .AND. FA(10).NE.0.0) THEN RECOUT(4) = (FA(5)+FA(9))*0.5 RECOUT(5) = SQRT(FA(6)**2+FA(10)**2)*0.5 RECOUT(6) = FA(5) - FA(9) RECOUT(7) = SQRT(FA(6)**2+FA(10)**2) C Case 2: I+ but not I- present. ELSEIF (FA(6).NE.0.0 .AND. FA(10).EQ.0.0) THEN ISYMB = 1 RECOUT(4) = FA(5) RECOUT(5) = FA(6) C Case 3: I- but not I+ present. ELSEIF (FA(6).EQ.0.0 .AND. FA(10).NE.0.0) THEN ISYMB = 2 RECOUT(4) = FA(9) RECOUT(5) = FA(10) C C---- We want to keep I+ and I- components if they exist. However, to C ensure backwards compatibility if I+ = I- = 0 when Imean > 0 set C I+ = I- = Imean. This should NEVER happen with new SCALA. C ELSE IF (.NOT.(LOGMSS(4) .OR. LOGMSS(5))) THEN IF (RECIN(5) .NE. 0.0) THEN RECOUT(4) = FA(1) RECOUT(5) = FA(2) RECOUT(8) = FA(1) RECOUT(9) = FA(2) RECOUT(14) = 0.01*RECIN(4)*SCALEF*SCALEF RECOUT(15) = 0.01*RECIN(5)*SCALEF*SCALEF RECOUT(10) = FA(1) RECOUT(11) = FA(2) RECOUT(16) = 0.01*RECIN(4)*SCALEF*SCALEF RECOUT(17) = 0.01*RECIN(5)*SCALEF*SCALEF END IF END IF END IF C if centric.... IF (ICEN.EQ.2) THEN RECOUT(6) = 0.0 RECOUT(7) = 0.0 ENDIF C RECOUT(18) = RECIN(16) RECOUT(NLPRGO) = REAL(ISYMB) END IF IF (.NOT.(QISNAN(RECOUT(4)).OR.QISNAN(RECOUT(5)))) THEN IF (RECOUT(4).LT.FMIN) THEN FMIN = RECOUT(4) SIGMIN = RECOUT(5) ELSE IF (RECOUT(4).GE.FMAX) THEN FMAX = RECOUT(4) SIGMAX = RECOUT(5) END IF END IF C C---- Find systematic absences C INHKL(1) = NINT(RECOUT(1)) INHKL(2) = NINT(RECOUT(2)) INHKL(3) = NINT(RECOUT(3)) C C ************************ CALL EPSLON(INHKL,WEIGHT,ISYSAB) C ************************ C C Store systematic absences: OMIT them from output IF (ISYSAB.GT.0) THEN NSYSABS = NSYSABS+1 IF (NSYSABS .LE. MAXSYSABS) THEN DO 222, JKD=1,3 IHKLSYSABS(JKD,NSYSABS) = INHKL(JKD) 222 CONTINUE IF (QISNAN(RECOUT(4)).OR.QISNAN(RECOUT(5))) THEN FSIGSYSABS(1,NSYSABS) = 0.0 FSIGSYSABS(2,NSYSABS) = 0.0 ELSE FSIGSYSABS(1,NSYSABS) = RECOUT(4) FSIGSYSABS(2,NSYSABS) = RECOUT(5) ENDIF GO TO 180 ENDIF ENDIF C C---- CALL to write a reflection record to an MTZ file which C has been opened for write. Also output column labels etc. C should already have been setup using LWASSN. C C ========================== CALL LWREFL(MTZOUT,RECOUT) C ========================== C IOUT = IOUT + 1 C End of output loop END IF C C GO TO 180 C C---- End of file reached on 2nd pass C 230 CONTINUE C C---- CALL to close an MTZ file which has been opened for read. C C ============= CALL LRCLOS(MTZIN) C ============= C C---- CALL to close an MTZ file which has been opened for write. C The new header information should already have been supplied by C calls to other LW* routines. This CALL writes the MTZ C header and its associated history header to the output file, C writes a pointer to the headers in the first record of the file C and closes the file. The headers are actually written at the end C of the file. C call ccp4h_pre_end() call ccp4h_rule() call ccp4h_header('Header Information to Output MTZ File', + 'header',3) call ccp4h_pre_beg() IF (LOGOUT) THEN CALL LWCLOS(MTZOUT,MTZPRT) ELSE WRITE(6,'(A)') + ' No output MTZ file written. ' ENDIF C call ccp4h_summary_beg() WRITE (LUNOUT,FMT=6020) NREF,IOUT 6020 FORMAT (////,' Number of reflections input = ',I10,/, + ' Number of terms output = ',I10,/) call ccp4h_summary_end() WRITE (LUNOUT,FMT=6021) NREJ 6021 FORMAT (/,'Number of terms rejected =',I10,/, + ' ( having EITHER Iobs .LT. -3.7*SDobs OR I', + 'obs .LT. (SDobs)**2/MeanI - 4.0*SDobs )') call ccp4h_pre_end() C C Print systematic absences C IF (NSYSABS .GT. 0) THEN WRITE(LUNOUT,FMT=6100) 6100 FORMAT(//,21X,'Systematic absences',/, + 21X,'===================',/) C WRITE(LUNOUT,'(11X,A,//,A,/)') + 'Systematic absences are OMITTED from output file', $ ' h k l F sd' DO 590, J=1,MIN(NSYSABS,MAXSYSABS) WRITE(LUNOUT,FMT=6110) $ (IHKLSYSABS(JJ,J),JJ=1,3),(FSIGSYSABS(JJ,J),JJ=1,2) 6110 FORMAT(1X,3I5,2F10.2) C C---- harvest C C loop_ C _refln_sys_abs.h C _refln_sys_abs.k C _refln_sys_abs.l C _refln_sys_abs.I C _refln_sys_abs.sigmaI C _refln_sys_abs.I_over_sigmaI C CALL Hrefln_sys_abs( IHKLSYSABS(1,J), + IHKLSYSABS(2,J), + IHKLSYSABS(3,J), + FSIGSYSABS(1,J), + FSIGSYSABS(2,J) ) C C---- harvest C 590 CONTINUE IF (NSYSABS .GT. MAXSYSABS) THEN WRITE(LUNOUT,'(1X,A,I6,A)') '***** ',NSYSABS-MAXSYSABS, $ ' systematic absences omitted from printed list *****' ENDIF ENDIF C C---- Output Wilson's distribution . call ccp4h_rule() call ccp4h_header('Distributions of Observed Intensity', + 'obsd',3) call ccp4h_pre_beg() C WRITE (LUNOUT,FMT=6022) 6022 FORMAT (/,' Distributions of Observed Intensity Magnitudes',/, + ' ----------------------------------------------',/) WRITE (LUNOUT,FMT=6024) 6024 FORMAT ( +' Tables below give percentage of terms for which I.le.Z ',/, +' where Z is defined as I/Mn(I) for the range of 4*((Sint', +'heta/Lamda)**2)',//, $ ' Also the 2nd, 3rd & 4th moments of I, Mn(I**k)/Mn(I)**k', $ ' for k = 2,3,4 (labelled Moment2, Moment3, Moment4))',// +' Z values in tables :',/, + 15X,'0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0') C C DO 250 JDO250 = 2,10 DO 240 JDO240 = 1,2 FRACT(JDO250,JDO240) = FRACT(JDO250-1,JDO240) + + FRACT(JDO250,JDO240) 240 CONTINUE 250 CONTINUE C C DO 280 JDO280 = 1,NRANGE DO 270 JDO270 = 2,10 DO 260 JDO260 = 1,2 FRACS(JDO280,JDO270,JDO260) = FRACS(JDO280,JDO270-1, + JDO260) + FRACS(JDO280,JDO270,JDO260) 260 CONTINUE 270 CONTINUE 280 CONTINUE call ccp4h_pre_end() call ccp4h_summary_beg() call ccp4h_header('Acentric Moments of Intensity', + 'acentric',3) call ccp4h_summary_end() call ccp4h_pre_beg() C C DO 320 JDO320 = 1,2 IF (JDO320.EQ.1) THEN call ccp4h_summary_beg() WRITE (LUNOUT,FMT=6026) 6026 FORMAT (/,' ACENTRIC MOMENTS OF INTENSITY',/, + ' ------------------------------ ',/) WRITE (LUNOUT,FMT=6998) THA, THMMTA 6998 FORMAT (' THEORETICAL Distribution ',/,13X,10F5.1, $ 5X,3F8.2,//, + ' Observed distribution in ranges of 4*((Sintheta/Lamda)**2)') call ccp4h_summary_end() call ccp4h_graph_beg(0,0) WRITE (LUNOUT,FMT=6030)SMAXWL,SMAXWL,SMAXWL ELSE IF (JDO320.EQ.2) THEN call ccp4h_pre_end() IF(NTOT(2) .NE.0) THEN call ccp4h_summary_beg() call ccp4h_header('Centric Moments of Intensity', + 'centric',3) call ccp4h_summary_end() call ccp4h_pre_beg() WRITE (LUNOUT,FMT=6028) 6028 FORMAT (/,' CENTRIC MOMENTS OF INTENSITY: 1Bar N(Z) ',/, + ' ----------------------------- ',/) WRITE (LUNOUT,FMT=6999) THC, THMMTC 6999 FORMAT (' THEORETICAL Distribution ',/,13X,10F5.1, $ 5X,3F8.2,//, + ' Observed distribution in ranges of 4*((Sintheta/Lamda)**2)') call ccp4h_graph_beg(0,0) WRITE (LUNOUT,FMT=6031) SMAXWL,SMAXWL,SMAXWL 6030 FORMAT (' Range Nref',23X,'N(Z) ',29X, $ 'Moment2 Moment3 Moment4',/,' ----- ----',23X, + '---- ',29X,3('------- '),//, + ' $TABLE:', + ' Acentric Moments of (I**k)/(I)**k for k = 0.5,1.5,2,3,4:',/, + ' $GRAPHS', + ' : 1st & 3rd Moment of E (Expected values = 0.866, 1.329', + 'Perfect twin = 0.94, 1.175):0|',f5.3,'x0|2:2,17,18:',/, + ' : 2nd Moment of I or z or E**2 ', + '(Expected value = 2, Perfect Twin = 1.5)', + ':0|',f5.3,'x0|5:2,14:',/, + ' : 3rd & 4th Moment of I or z (Expected value = 6, 24', + ' Perfect Twin 3, 7.5):0|',f5.3,'x0|48:2,15,16:',/, + ' $$',/, + ' Resln_Range 1/resol^2 Nref Nz1 Nz2 Nz3 Nz4 Nz5 Nz6 ', + 'Nz7 Nz8 Nz9 Nz10', + ' MomentZ2 MomentZ3 MomentZ4 MomentE1 MomentE3 $$',/, '$$') 6031 FORMAT (' Range Nref',23X,'N(Z) ',29X, + 'Moment2 Moment3 Moment4',/,' ----- ----',23X, + '---- ',29X,3('------- '),//, + ' $TABLE:', + ' Centric Moments of (I**k)/(I)**k for k = 0.5,1.5,2,3,4:',/, + ' $GRAPHS', + ' : 1st & 3rd Moments (E) (Expected = 0.798,1.596)', + ':0|',f5.3,'x0|4:2,17,18:',/, + ' : 2nd Moment of I or z or E**2 (Expected = 3)', + ':0|',f5.3,'x0|5:2,14:',/, + ' : 3rd & 4th Moment of I or z (Expected = 15,105)', + ':0|',f5.3,'x0|120:2,15,16:',/, + ' $$',/, + ' Resln_Range 1/resol^2 Nref Nz1 Nz2 Nz3 Nz4 Nz5 Nz6 ', + 'Nz7 Nz8 Nz9 Nz10', + ' MomentZ2 MomentZ3 MomentZ4 MomentE1 MomentE3 $$',/, '$$') ENDIF ELSE GO TO 320 END IF C NMMTNM = 0 DO 321, K = 1,MAXMMT SMMTOT(K) = 0.0 SMMTOTE(K) = 0.0 321 CONTINUE C C DO 300 JDO300 = 1,NRANGE IF (NW(JDO300,JDO320).NE.0) THEN C NTOT(JDO320) = NTOT(JDO320) + NW(JDO300,JDO320) FN = NW(JDO300,JDO320) C C DO 290 JDO290 = 1,10 FRACS(JDO300,JDO290,JDO320) = (FRACS(JDO300,JDO290, + JDO320)/FN)*100.0 290 CONTINUE C IF (NMMNUM(JDO300,JDO320).GT.0 .AND. + SMMIMM(1,JDO300,JDO320).GT.0.0) THEN NMMTNM = NMMTNM + NMMNUM(JDO300,JDO320) DO 291, K=1,MAXMMT C for this resolution shell SMMIMM(K,JDO300,JDO320) = $ SMMIMM(K,JDO300,JDO320)/NMMNUM(JDO300,JDO320) SMMEMM(K,JDO300,JDO320) = $ SMMEMM(K,JDO300,JDO320)/NMMNUM(JDO300,JDO320) 291 CONTINUE C SIMEAN = SMMIMM(1,JDO300,JDO320) DO 292, K=1,MAXMMT C /**K for this resolution shell IF(K.NE.1)SMMIMM(K,JDO300,JDO320) = $ SMMIMM(K,JDO300,JDO320)/(SIMEAN**K) SMMTOT(K) = SMMTOT(K) + $ SMMIMM(K,JDO300,JDO320)*NMMNUM(JDO300,JDO320) IF(K.EQ.1 .OR. K.EQ.3) THEN SMMEMM(K,JDO300,JDO320) = $ SMMEMM(K,JDO300,JDO320)/(SIMEAN**(0.5*K)) SMMTOTE(K) = SMMTOTE(K) + $ SMMEMM(K,JDO300,JDO320)*NMMNUM(JDO300,JDO320) END IF 292 CONTINUE ELSE DO 293, K=2,MAXMMT SMMIMM(K,JDO300,JDO320) = 0.0 SMMEMM(K,JDO300,JDO320) = 0.0 293 CONTINUE ENDIF SS = (REAL(JDO300)*RANGE+SMINWL) WRITE (LUNOUT,FMT=6032) JDO300,SS,NW(JDO300,JDO320), + (FRACS(JDO300,JLOP,JDO320),JLOP=1,10), + (SMMIMM(K,JDO300,JDO320),K=2,MAXMMT), + SMMEMM(1,JDO300,JDO320),SMMEMM(3,JDO300,JDO320) 6032 FORMAT (1X,I5, F8.4,I7,10F7.1,2X,5F8.2) END IF 300 CONTINUE WRITE (LUNOUT,'(A)') ' $$ ' C call ccp4h_graph_end() C DO 310 JDO310 = 1,10 IF (NTOT(JDO320).GT.0) FRACT(JDO310, + JDO320) = (FRACT(JDO310,JDO320)/REAL(NTOT(JDO320)))*100.0 310 CONTINUE C DO 311, K=1,MAXMMT C </**K> for all resolutions IF (NMMTNM .GT. 0) THEN SMMTOT(K) = SMMTOT(K)/NMMTNM SMMTOTE(K) = SMMTOT(K)/NMMTNM ELSE SMMTOT(K) = 0.0 SMMTOTE(K) = 0.0 ENDIF 311 CONTINUE C call ccp4h_summary_beg() IF (NTOT(JDO320).GT.0) WRITE (LUNOUT,FMT=6034) NTOT(JDO320), + (FRACT(JLOP,JDO320),JLOP=1,10),(SMMTOT(K),K=2,MAXMMT), + SMMTOTE(1),SMMTOTE(3) 6034 FORMAT (/,' Totals of Observed Distributions (averages) :',/, $ 12X,I9,10F7.1,2X,5F8.2,/) C Is the second moment nearer 1.5 (twinning) than 2 (normal) C or the third moment nearer 3.0 (twinning) than 6 (normal) IF(KYTWIN .AND. JDO320.EQ.1 .AND. $ (ABS(SMMTOT(2)/1.5 -1.0) .LT. 0.2) .AND. $ (ABS(SMMTOT(3)/3.0 -1.0) .LT. 0.3) ) $ CALL CCPERR(2,' **** Beware! - Possible Twinning ****') C This needs a CCP_ERR message.. call ccp4h_summary_end() 320 CONTINUE call ccp4h_pre_end() C C---- Write distribution as table for plotting C call ccp4h_rule() call ccp4h_header('Cumulative Intensity Distribution', + 'cumm',3) call ccp4h_pre_beg() call ccp4h_graph_beg(0,0) WRITE (LUNOUT,FMT=6038) 6038 FORMAT(/, $ ' $TABLE: Cumulative intensity distribution:',/, $ ' $GRAPHS',/, $ ' :Cumulative intensity distribution (Acentric and centric)',/, $ ' :N:1,2,3,4,5:',/, $ ' $$',/ $ ' Z N(Z)Atheor N(Z)Acen N(Z)Ctheor N(Z)Cen $$', $/,' $$',/) ZERO = 0.0 WRITE(LUNOUT,FMT=6039) ZERO,ZERO,ZERO,ZERO,ZERO DO 325 I=1,10 WRITE (LUNOUT,FMT=6039) $ REAL(I)*0.1,THA(I),FRACT(I,1),THC(I),FRACT(I,2) 6039 FORMAT(1X,F10.1,4F10.1) C C---- harvest C C loop_ C _reflns_intensity_shell.Z C _reflns_intensity_shell.NZ_acentric_theory C _reflns_intensity_shell.NZ_acentric_observed C _reflns_intensity_shell.NZ_centric_theory C _reflns_intensity_shell.NZ_centric_observed C CALL Hreflns_intensity_shell(REAL(I)*0.1,THA(I), + FRACT(I,1),THC(I),FRACT(I,2)) C C---- harvest C 325 CONTINUE WRITE(LUNOUT,FMT='(/,A,//)') ' $$' call ccp4h_graph_end() call ccp4h_pre_end() call ccp4h_rule() C C C---- Print sd statistics TOTIN = 0.0 TOTII = 0.0 TOTIS = 0.0 TOTIR = 0.0 TOTFI = 0.0 TOTFS = 0.0 TOTFR = 0.0 C call ccp4h_header('Mean Amplitude vs. Resolution', + 'resolution',3) call ccp4h_pre_beg() call ccp4h_graph_beg(0,0) WRITE(LUNOUT,6400) 6400 FORMAT(/, $' $TABLE: Amplitude analysis against resolution:',/, $' $GRAPHS:Mn(F) v. resolution:N:2,9::', $'Mn(F/sd) v. resolution:N:2,12:',/,' $$',/, $ ' Range 1/resol^2 Dmax Nref Mn(I) Mn(sd) Mn(I)/Mn(sd) ', $ ' Mn(I/sd) Mn(F) Mn(sd) Mn(F)/Mn(sd) Mn(F/sd) $$',/,' $$') C DO 400, I=1,KNBINS IF (SUMINS(I) .GT. 0.000001) THEN TOTIN = TOTIN + SUMINS(I) TOTII = TOTII + SUMIIS(I) TOTIS = TOTIS + SUMISS(I) TOTIR = TOTIR + SUMIRS(I) TOTFI = TOTFI + SUMFIS(I) TOTFS = TOTFS + SUMFSS(I) TOTFR = TOTFR + SUMFRS(I) SS = (REAL(I)*RANGE+SMINWL) DS = 1./SQRT(SS) XJ=SUMIRS(I)/SUMINS(I) AIMEAN=SUMIIS(I)/SUMINS(I) SIGIMEAN=SUMISS(I)/SUMINS(I) XF=SUMFRS(I)/SUMINS(I) AFMEAN=SUMFIS(I)/SUMINS(I) SIGFMEAN=SUMFSS(I)/SUMINS(I) J=SUMINS(I) WRITE(LUNOUT,6410) I,SS,DS,J, $ AIMEAN,SIGIMEAN,SUMIIS(I)/SUMISS(I),XJ, $ AFMEAN,SIGFMEAN,SUMFIS(I)/SUMFSS(I),XF 6410 FORMAT(1X,I4,2X,F10.4,2X,F10.2,2X,I8,2(4F12.2,2X)) C C---- harvest C R2 = DS IF (I.eq. 1) THEN R1 = RMINWL ELSE R1 = (REAL(I-1)*RANGE+SMINWL) R1 = 1./SQRT(R1) END IF C C loop_ C _reflns_scaling_shell.d_res_high C _reflns_scaling_shell.d_res_low C _reflns_scaling_shell.num_reflns_observed C _reflns_scaling_shell.mean_obs C _reflns_scaling_shell.mean_obs C CALL Hreflns_scaling_shell(R1,R2,J,XJ,XF) C C---- harvest C ENDIF 400 CONTINUE C WRITE(LUNOUT,'(A)') ' $$' call ccp4h_graph_end() C XJ=TOTIR/TOTIN AIMEAN=TOTII/TOTIN SIGIMEAN=TOTIS/TOTIN XF=TOTFR/TOTIN AFMEAN=TOTFI/TOTIN SIGFMEAN=TOTFS/TOTIN J=TOTIN WRITE(LUNOUT,6420) J, $ AIMEAN,SIGIMEAN,TOTII/TOTIS,XJ, $ AFMEAN,SIGFMEAN,TOTFI/TOTFS,XF 6420 FORMAT(/,1X,' TOTALS ',I8,2(4F9.2)) C C---- harvest C C _reflns.B_iso_Wilson_estimate WILBT(1) C _reflns.entry_id Projectname C _reflns.dataset_id Lead_acetate_derivative C _reflns.d_resolution_high RMINWL C _reflns.d_resolution_low RMAXWL C _reflns.observed_criterion C 'reject either I<-3.7*sigmaI or I<(sigmaI)**2/MeanI-4.0*sigmaI' C _reflns.mean_obs_all XJ C _reflns.mean_obs_all XF C _reflns.number_obs J C CALL Hreflns(WILBT(1),RMINWL,RMAXWL, +'reject either I<-3.7*sigmaI or I<(sigmaI)**2/MeanI-4.0*sigmaI', + XJ,XF,J) C C---- harvest C C FMIN etc. generated from RECOUT IF(LOGOUT) WRITE (LUNOUT,FMT=6036) FMIN,SIGMIN,FMAX,SIGMAX 6036 FORMAT (///,' Minimum F = ',F10.3,/, + ' with SD = ',F10.3,/, + ' Maximum F = ',F10.3,/, + ' with SD = ',F10.3,//) call ccp4h_pre_end() call ccp4h_rule() C C----Display anisotropy with FALLOFF module from Yorgo Modis. call ccp4h_header('Anisotropic Analysis: FALLOFF','anis',3) call ccp4h_pre_beg() IF (KYFALL) THEN CALL FALLOFF ENDIF C c-harvest C C----- every program must call Hclose C only with this call is the deposit file written C CALL Hclose C C---harvest C CALL CCPERR(0,'Normal termination') END C C FUNCTION GAMMLN(XX) C =================== C C---- Returns the value ln(G(xx) for xx > 0 . Full accuracy for C xx>1.0 C C Use formula G(1-z) = PI*z/( G(1+z)*sin(pi*z) ) C to get full accuracy. C C IMPLICIT NONE C INTEGER J REAL XX,GAMMLN DOUBLE PRECISION COF(6),STP,HALF,ONE,FPF,X,TMP,SER C C DATA COF,STP/76.18009172D0,-86.50532033D0,24.01409822D0, 1 -1.231739516D0,0.120858003D-2,-0.536382D-5, 1 2.50662827465D0/ C DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/ C C X=XX-ONE TMP=X+FPF TMP=(X+HALF)*LOG(TMP)-TMP SER=ONE DO 11 J=1,6 X=X+ONE SER=SER+COF(J)/X 11 CONTINUE GAMMLN=TMP+LOG(STP*SER) RETURN END C C C FUNCTION GAMMP(A,X) C =================== C C C---- Returns the incomplete GAMMA function Q(a,x) = 1-P(a,x) C REAL GAMMP,GAMMCF,GLN,GAMSER,A,X C IF (X.LT.0.0 .AND. A.LE.0.0) $ CALL CCPERR(1,' Garbage - x,a lt 0 ') IF(X.LT.A+1.0) THEN CALL GSER(GAMSER,A,X,GLN) GAMMP = GAMSER ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMP = 1.0 - GAMMCF ENDIF RETURN END C C C FUNCTION GAMMQ(A,X) C =================== C C---- Returns the incomplete GAMMA function Q(a,x) = 1-P(a,x) C REAL GAMMQ,GAMMCF,GLN,GAMSER,A,X C IF (X.LT.0.0 .AND. A.LE.0.0) $ CALL CCPERR(1,' Garbage - x,a lt 0 ') IF(X.LT.A+1.0) THEN C CALL GSER(GAMSER,A,X,GLN) C GAMMQ = 1.0 - GAMSER ELSE C CALL GCF(GAMMCF,A,X,GLN) C GAMMQ = GAMMCF ENDIF RETURN END C C C SUBROUTINE GCF(GAMMCF,A,X,GLN) C =============================== C C---- Returns the incomplete GAMMA function Q(a,x) C evaluated by its continued fraction representation as GAMMCF C INTEGER ITMAX REAL EPS PARAMETER (ITMAX=100,EPS=3.E-7) C REAL GLN,GAMMLN,A,GOLD,A0,A1,X,B0,B1,FAC,AN,ANA,ANF,G,GAMMCF INTEGER N C GLN = GAMMLN(A) GOLD = 0.0 A0= 1.0 A1= X B0= 0.0 B1= 1.0 FAC= 1.0 C C DO 11 N=1,ITMAX AN = REAL(N) ANA = AN - A A0=(A1+A0*ANA)*FAC B0=(B1+B0*ANA)*FAC ANF = AN*FAC A1=X*A0 + ANF*A1 B1=X*B0 + ANF*B1 C C IF(A1.NE.0.0) THEN FAC = 1.0/A1 G = B1*FAC IF( ABS((G-GOLD)/G).LT.EPS) GO TO 1 GOLD = G ENDIF 11 CONTINUE C C CALL CCPERR(1,' ** A too large - itmax too small***') C 1 GAMMCF = EXP(-X + A*ALOG(X) -GLN)*G RETURN END C C C SUBROUTINE GSER(GAMSER,A,X,GLN) C ============================= C C C---- Returns the incomplete GAMMA function 1-P(a,x) C evaluated by its series representation as GAMSER. C C INTEGER ITMAX REAL EPS PARAMETER (ITMAX=100,EPS=3.E-7) C REAL GLN,GAMMLN,A,X,GAMSER,AP,SUM,DEL INTEGER N C C GLN = GAMMLN(A) C C IF (X.LE.0.0)THEN IF (X.LT.0.0) CALL CCPERR(1, ' Garbage - x lt 0 ') GAMSER = 0.0 RETURN ENDIF C C AP = A SUM = 1.0/A DEL = SUM C C DO 11 N=1,ITMAX AP = AP + 1.0 DEL = DEL*X/AP SUM = SUM + DEL IF(ABS(DEL).LT.ABS(SUM)*EPS) GO TO 1 11 CONTINUE C C CALL CCPERR(1,' ** A too large - itmax too small***') C 1 GAMSER = SUM*EXP(-X + A*LOG(X) -GLN) RETURN END C C C SUBROUTINE ITRANS(KINT,KCHAR) C ============================ C C C---- Translate integers into character strings. C C .. Scalar Arguments .. INTEGER KINT C .. C .. Array Arguments .. CHARACTER KCHAR(5)*1 C .. C .. Scalars in Common .. REAL CONV INTEGER ISYSR,ISYSW,NCOL,NMON 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 .. Common blocks .. COMMON /IO/ISYSR,ISYSW,NCOL,CONV,NMON C .. C .. Data statements .. DATA BLANK/' '/,MINUS/'-'/ DATA ITABLE/'0','1','2','3','4','5','6','7','8','9'/ C .. C C NEG = .FALSE. IF (KINT.LT.0) NEG = .TRUE. C C---- Initialize KCHAR to represent 0, in case it is. C DO 10 I = 1,4 KCHAR(I) = BLANK 10 CONTINUE KCHAR(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 KCHAR(J) = ITABLE(I) J = J - 1 INT2 = INT(INT2/10.0) GO TO 20 END IF IF (NEG) KCHAR(J) = MINUS END C C C SUBROUTINE HKLLIM + (NLAUE,DMIN,DMAX,CELL,HMIN,HMAX,KMIN,KMAX,LMIN,LMAX) C =============================================================== C C C 22/3/93 Changed to match Laue groups from ASUSET, but code left c here for some other numbers which are not used by ASUSET, c viz 1,2,5,16 Phil Evans C C C .. Scalar Arguments .. REAL DMAX,DMIN INTEGER HMAX,HMIN,KMAX,KMIN,LMAX,LMIN,NLAUE C .. C .. Arrays in Common .. C REAL CELL(6),RABC,RABCAN,RABCSQ,RANGD REAL CELL(6) C .. C .. Intrinsic Functions .. INTRINSIC INT C .. C .. Save statement .. SAVE C .. C C IF (NLAUE.LT.1 .OR. NLAUE.GT.16) GO TO 10 C C---- otherwise set them from resolution command C C Note this IS CORRECT ! C HMAX = INT(CELL(1)/DMAX) KMAX = INT(CELL(2)/DMAX) LMAX = INT(CELL(3)/DMAX) C HMIN = 0 KMIN = 0 LMIN = 0 C C IF (NLAUE.EQ.1) THEN C C---- LAUE = 1, 1bar, hkl:h>=0 0kl:k>=0 00l:l>=0 C*C Not used by ASUSET C KMIN = -KMAX LMIN = -LMAX RETURN C C---- LAUE = 2, 1bar, hkl:k>=0 h0l:l>=0 h00:h>=0 C*C Not used by ASUSET C ELSE IF (NLAUE.EQ.2) THEN HMIN = -HMAX LMIN = -LMAX RETURN C C---- LAUE = 3, 1bar, hkl:l>=0 hk0:h>=0 0k0:k>=0 C ELSE IF (NLAUE.EQ.3) THEN HMIN = -HMAX KMIN = -KMAX RETURN C C---- LAUE = 4, 2/m, hkl:k>=0, l>=0 hk0:h>=0 C ELSE IF (NLAUE.EQ.4) THEN HMIN = -HMAX RETURN C C---- LAUE = 5, 2/m, hkl:h>=0, l>=0 0kl:k>=0 (2-nd sett.) C*C Not used by ASUSET C ELSE IF (NLAUE.EQ.5) THEN KMIN = -KMAX RETURN C C---- LAUE = 6, mmm, hkl:h>=0, k>=0, l>=0 C ELSE IF (NLAUE.EQ.6) THEN RETURN C C---- LAUE = 7, 4/m, hkl:h>=0, k>0, l>=0 with k>=0 for h=0 C ELSE IF (NLAUE.EQ.7) THEN RETURN C C---- LAUE = 8, 4/mmm, hkl:h>=0, h>=k>=0, l>=0 C ELSE IF (NLAUE.EQ.8) THEN RETURN C C---- LAUE = 9, 3bar, hkl:h>=0, k>0 including 00l:l>0 C ELSE IF (NLAUE.EQ.9) THEN KMAX = HMAX LMIN = -LMAX C Cejd - Loses all 00l C KMIN = 1 RETURN C C---- LAUE = 10, 312, hkl:h>=0, k>=0 with k<=h and l>=0 if h = 0 C ELSE IF (NLAUE.EQ.10) THEN KMAX = HMAX LMIN = -LMAX RETURN C C---- LAUE = 11, 321, hkl:h>=0, k>=0 with k<=h and l>=0 if h = k C ELSE IF (NLAUE.EQ.11) THEN KMAX = HMAX LMIN = -LMAX RETURN C C---- LAUE = 12, 6/m, hkl:h>=0, k>0, l>=0 with k>=0 for h=0 C ELSE IF (NLAUE.EQ.12) THEN KMAX = HMAX RETURN C C---- LAUE = 13, 6/mmm, hkl:h>=0, h>=k>=0, l>=0 C ELSE IF (NLAUE.EQ.13) THEN KMAX = HMAX RETURN C C---- LAUE = 14, m3, hkl:h>=0, k>=0, l>=0 with l>=h, k>=h for l=h C and k>h for l>h C ELSE IF (NLAUE.EQ.14) THEN RETURN C C---- LAUE = 15, m3m, hkl:k>=l>=h>=0 C ELSE IF (NLAUE.EQ.15) THEN RETURN C C---- LAUE = 16, m3m, hkl:h>=k>=l>=0 C*C Not used by ASUSET C ELSE IF (NLAUE.EQ.16) THEN RETURN END IF C C 10 WRITE (6,FMT=6000) NLAUE 6000 FORMAT (//,1X,'**** WRONG LAUE GROUP = ',I6,' ****',/) C C ****************************** CALL CCPERR(1,' wrong laue group ') C ****************************** C C END C C SUBROUTINE KEYIN C ================ C C IMPLICIT NONE C C C C .. Parameters .. INTEGER NPAR PARAMETER (NPAR=200) INTEGER MAXSER PARAMETER (MAXSER=60) INTEGER MAXSYM PARAMETER (MAXSYM=96) INTEGER NPARM PARAMETER (NPARM=200) INTEGER NLPRGI PARAMETER (NLPRGI=16) INTEGER ncols PARAMETER (ncols=19) INTEGER MAXHIS PARAMETER (MAXHIS=30) C .. C .. Scalars in Common .. REAL RANGE,cone INTEGER ISOUT,KANOM,LABFLG,LBOTKN,LENTIT,LUNIN,LUNOUT,MTITOT, + MTZBPR,MTZERR,MTZFLG,MTZIN,MTZOUT,MTZPRT,NBATX, + NHSOUT,NLAUE,NREFLX,NSPGRP,NSYM,NSYMP, + NUMCOL,NRANGE LOGICAL CSYMM,KYCELL,KYHIST,KYSYMM,KYTITL,KYRESO,LBOLOG, $ LNOTRN,MTZEOF,RNGSET,KYFALL,PLOTX,PLOTFALL CHARACTER LATTYP*1,DATEMT*8,LAUNAM*10,NAMSPG*10,PGNAME*10, + VERSNX*10,TITIN*70,TITOUT*70,LBOLIN*400 C .. C .. Arrays in Common .. REAL CELL,CELMTZ,CELOUT,RNGMTZ,RSYM,RSYMT INTEGER LBOBEG,LBOEND,MTZLOK CHARACTER CTPRGO*1,CTPRGI*1,LSPRGO*30, + LSPRGI*30,HISOUT*80 C .. C .. Local Scalars .. INTEGER IPRINT,J,JDO,LCOUNT,LXX,NTOK,I, $ IAT,NRES LOGICAL LEND CHARACTER KEY*4,CWORK*400,LINE*400 C .. C .. Local Arrays .. REAL FVALUE(NPAR) INTEGER IBEG(NPAR),IDEC(NPAR),IEND(NPAR),ITYP(NPAR) CHARACTER CVALUE(NPAR)*4 C .. C .. External Functions .. INTEGER LENSTR EXTERNAL LENSTR C .. C .. External Subroutines .. EXTERNAL CCPUPC,KEYNUM,PARSER,RDSYMM C .. C .. Intrinsic Functions .. INTRINSIC LEN,NINT C .. C .. Common blocks .. COMMON /AGRSYM/NSPGRP,NSYM,NSYMP,NLAUE,CSYMM, + RSYM(4,4,MAXSYM),RSYMT(4,4,MAXSYM) COMMON /CHRMTZ/LSPRGI(16),CTPRGI(16), + LSPRGO(ncols),CTPRGO(ncols),HISOUT(MAXHIS), + VERSNX,TITIN,TITOUT,DATEMT,LBOLIN COMMON /CHRSYM/PGNAME,NAMSPG,LAUNAM,LATTYP COMMON /INOUT/LUNIN,LUNOUT,ISOUT COMMON /KEYISC/KANOM,LNOTRN COMMON /KEYRAR/CELL(6) COMMON /KEYRSC/RANGE,NRANGE COMMON /WRKMTZ/MTZIN,MTZOUT,MTZERR,MTZBPR,MTZPRT,MTZFLG,MTITOT, + NBATX,NHSOUT,NUMCOL,NREFLX,LENTIT,LABFLG,MTZEOF, + KYTITL,KYHIST,KYCELL,KYSYMM,KYRESO,MTZLOK(NLPRGI), + CELMTZ(6),CELOUT(6),RNGMTZ(2,16), + LBOBEG(NPARM),LBOEND(NPARM),LBOLOG,LBOTKN COMMON /RNG/RNGSET COMMON /FALL/KYFALL,PLOTX,PLOTFALL,cone C COMMON /APLSCL/ SCALEF, LAPSCL REAL SCALEF LOGICAL LAPSCL C .. C-harvest CHARACTER PNAME_KEYWRD*64,DNAME_KEYWRD*64 COMMON /CHARV/ PNAME_KEYWRD,DNAME_KEYWRD LOGICAL PRIVATE,USECWD,NOHARVEST LOGICAL GOT_PNAME,GOT_DNAME INTEGER NROWLIMIT,IvalueNotDet REAL ValueNotDet COMMON /IHARV/PRIVATE,USECWD,NOHARVEST, + GOT_PNAME,GOT_DNAME, + NROWLIMIT,IvalueNotDet, + ValueNotDet C-harvest C C wilson bits C INTEGER NBIN,NMON,IPLFSQ,NK REAL RMAXWL,RMINWL,SMINWL,SMAXWL,SMINSC,SMAXSC,RMINSC,RMAXSC, $ VPA LOGICAL SCBLOG,LOGCON,LOGRES,KYTWIN INTEGER IWITCH COMMON /WILKEY/ NBIN,NMON,RMAXWL,RMINWL,SMINWL,SMAXWL,SMINSC, + SMAXSC,RMINSC,RMAXSC,VPA, + IPLFSQ,NK COMMON /LOGWIL/ IWITCH,SCBLOG,KYTWIN C CHARACTER NW*2 CHARACTER NW*4 COMMON /WILCHR/ NW(20) INTEGER NC,KatmWghts REAL FormFactAcoeff,FormFactBcoeff COMMON /WILFRM/ NC(20),KatmWghts(20),FormFactAcoeff(5,20), + FormFactBcoeff(4,20) C C C .. C C Default apply Wilson scale LAPSCL = .TRUE. SCALEF = 0.0 C GOT_PNAME = .false. GOT_DNAME = .false. NHSOUT = 0 NK = 0 LBOLOG = .FALSE. KYSYMM = .FALSE. KYTITL = .FALSE. KYRESO = .FALSE. SCBLOG = .FALSE. IWITCH = 1 LOGCON = .FALSE. LOGRES = .FALSE. IPRINT = 0 LCOUNT = 0 KANOM = -1 NRANGE = 60 RANGE = -1.0 NBIN = 60 VPA = 10.0 IPLFSQ = 0 NMON = 1000 CELL(1) = 0.0 KYCELL = .false. KYFALL = .TRUE. CONE = 30.0 plotx=.true. plotfall=.false. C c-harvest C Init variables for deposit files C PRIVATE = .false. USECWD = .false. NOHARVEST = .false. PNAME_KEYWRD = ' ' DNAME_KEYWRD = ' ' NROWLIMIT = 80 C C 10 CONTINUE LINE = ' ' KEY = ' ' 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 50 C C C---- Get chars of 1st token C KEY = LINE(IBEG(1) :IEND(1)) KEY = KEY(1:4) C C---- Comment !(Extra check also in Parse) C IF (KEY(1:1).EQ.'!') GO TO 10 IF (KEY(1:3).EQ.'END') GO TO 50 C C---- Convert to Uppercase C C *********** CALL CCPUPC(KEY) C *********** C C DATA (LSPRGI(J),J=1,13)/'H','K','L','IMEAN', C + 'SIGIMEAN','IW(+)','SIGIW(+)','I(+)','SIGI(+)', C + 'IW(-)','SIGIW(-)','I(-)','SIGI(-)'/ C---- LABI [optional input] C ========================== C IF (KEY.EQ.'LABI') THEN C ***************************************** CALL LKYIN(1,LSPRGI,NLPRGI,NTOK,LINE,IBEG,IEND) C ***************************************** C Range is the width of ranges on 4sin**2theta/lambda**2 C for the intensity analysis, OR the number of ranges GO TO 10 C C C---- RANGE n [optional input] C ========================== C C Range is the width of ranges on 4sin**2theta/lambda**2 C for the intensity analysis, OR the number of ranges C ELSE IF (KEY.EQ.'RANG') THEN C C IF (NTOK.LE.1) THEN NRANGE = 60 WRITE (LUNOUT,FMT=6000) 6000 FORMAT (//,' **** KEYWORDING WARNING ****',/, + ' No argument given on RANGE command,', $ ' using default = 60 bins') ELSE C C ************************************ CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C RANGE = FVALUE(2) RNGSET=.TRUE. IF (RANGE .GT. 1.0) THEN NRANGE = MIN(60,NINT(FVALUE(2))) C Set RANGE .lt. 0 as flag that NRANGE has been set RANGE = -1.0 ENDIF C END IF GO TO 10 C C---- CELL unit cell dimensions 3 or 6 real values [optional input] C =============================================================== C C if used overides cell dimensions in MTZ file C a, b, c, alpha, beta, gamma C C Cell dimensions in Angstroms and degrees C (angles default to 90.0). C If not given then the cell dimensions are taken from C the input MTZ file. C ELSE IF (KEY.EQ.'CELL') THEN CALL RDCELL(2,ITYP,FVALUE,NTOK,CELL) WRITE (LUNOUT,FMT=6004) 6004 FORMAT (//,' **** CELL Dimensions in MTZ file ignored ****') KYCELL = .true. GO TO 10 C C---- TRUNCATE "string" [OPTIONAL INPUT] default YES C ================================================= C C Use one of following C TRUNCATE YES (default) C or C TRUNCATE NO C ELSE IF (KEY.EQ.'TRUN') THEN C C IF (NTOK.LE.1) THEN WRITE (LUNOUT,FMT=6006) 6006 FORMAT (//,' **** KEYWORD WARNING ****',/, + ' No argument given on TRUNCATE KEY using default of YES') LNOTRN = .FALSE. ELSE IF (ITYP(2).EQ.2) THEN WRITE (LUNOUT,FMT=6008) 6008 FORMAT (//,' **** KEYWORD WARNING ****',/, + ' Integer argument given on TRUNCATE KEY using default of YES') LNOTRN = .FALSE. ELSE CWORK = LINE(IBEG(2) :IEND(2)) C C ************* CALL CCPUPC(CWORK) C ************* C IF (CWORK(1:3).EQ.'YES') LNOTRN = .FALSE. IF (CWORK(1:2).EQ.'NO') LNOTRN = .TRUE. END IF GO TO 10 C C---- ANOMALOUS "string" [OPTIONAL INPUT] default NO C ================================================= C C Use one of following C ANOMALOUS NO (default) C or C ANOMALOUS YES C ELSE IF (KEY.EQ.'ANOM') THEN C C IF (NTOK.LE.1) THEN WRITE (LUNOUT,FMT=6010) 6010 FORMAT (//,' **** KEYWORD WARNING ****',/, + ' No argument given on ANOMALOUS KEY using default of NO') KANOM = 0 ELSE IF (ITYP(2).EQ.2) THEN WRITE (LUNOUT,FMT=6012) 6012 FORMAT (//,' **** KEYWORD WARNING ****',/, + ' Integer argument given on ANOMALOUS KEY using default of NO') KANOM = 0 ELSE CWORK = LINE(IBEG(2) :IEND(2)) C C ************* CALL CCPUPC(CWORK) C ************* C IF (CWORK(1:3).EQ.'YES') KANOM = 1 IF (CWORK(1:2).EQ.'NO') KANOM = 0 END IF GO TO 10 C C C---- TITLE "string" [OPTIONAL INPUT] C ================================ C C C C Title defaults to 'FROM TRUNCATE' C C otherwise the input string is used C Title (max 80 characters) for the output MTZ file C ELSE IF (KEY.EQ.'TITL') THEN IF (NTOK.LE.1) GOTO 10 LXX = LEN(LINE) IF (LXX.GT.80) LXX = 80 TITOUT = LINE(IBEG(2) :LXX) KYTITL = .TRUE. GO TO 10 c c---- FALL YES|NO CONE PLTX|PLTY [OPTIONAL INPUT for FALLOFF] c ================================= c else if (key.eq.'FALL') then IF (NTOK.LE.1) THEN WRITE (LUNOUT,FMT=6106) 6106 FORMAT (//,' **** KEYWORD WARNING ****',/, + ' No argument given on FALLOFF KEY: using defaults') GO TO 10 ELSE CWORK = LINE(IBEG(2) :IEND(2)) CALL CCPUPC(CWORK) IF (CWORK(1:3).EQ.'YES') THEN KYFALL = .TRUE. ELSEIF (CWORK(1:2).EQ.'NO') THEN KYFALL = .FALSE. ELSE WRITE (LUNOUT,FMT=6108) 6108 FORMAT (//,' **** KEYWORD WARNING ****',/, +' Second argument of FALLOFF must be YES or NO: using defaults') GO TO 10 END IF END IF N = 3 6110 IF (N.GT.NTOK) GO TO 10 CWORK = LINE(IBEG(N) :IEND(N)) CALL CCPUPC(CWORK) IF (CWORK(1:4).EQ.'CONE') THEN N = N + 1 IF (N.LE.NTOK) CALL GTPREA(N,CONE,NTOK,ITYP,FVALUE) N = N + 1 GO TO 6110 ELSEIF (CWORK(1:4).EQ.'PLTX') THEN plotx=.TRUE. plotfall=.true. N = N + 1 GO TO 6110 ELSEIF (CWORK(1:4).EQ.'PLTY') THEN plotx=.FALSE. plotfall=.true. N = N + 1 GO TO 6110 ELSE WRITE (LUNOUT,FMT=6112) 6112 FORMAT (//,' **** KEYWORD WARNING ****',/, + ' Unrecognised subkeyword of FALLOFF KEY: ignoring') N = N + 1 GO TO 6110 END IF GO TO 10 C C---- HEADER keywords [OPTIONAL INPUT] C ================================= C Keywords to set flags for printing MTZ header C Subkeywords: NONE, BRIEF, HISTORY, ALL, NOBATCH, BATCH, ORIENTATION ELSE IF (KEY.EQ.'HEAD') THEN CALL RDHEAD(2,LINE,IBEG,IEND,ITYP,FVALUE,NTOK, . MTZPRT,MTZBPR) GO TO 10 C C---- LABOUT keyword C ============== C C C---- LABOUT in form F=userlabel etc [OPTIONAL INPUT] C ================================================ C C The labels KEY is optional C C Labels allowed are F SIGF DANO SIGDANO ISYM C C If this KEY is blank or absent, the default labels C F SIGF DANO SIGDANO ISYM are used. C C Note that any or all labels may be changed from default values C C C---- save it C ELSE IF (KEY.EQ.'LABO') THEN IF (.NOT.LBOLOG) THEN LBOLOG = .TRUE. LBOLIN = LINE LBOTKN = NTOK C C DO 40 JDO = 1,LBOTKN LBOBEG(JDO) = IBEG(JDO) LBOEND(JDO) = IEND(JDO) 40 CONTINUE C C GO TO 10 ELSE C C WRITE (LUNOUT,FMT=*) + ' **** ERROR Sorry only one LABOUT line allowed' CALL CCPERR(1,'Superfluous LABOUT command') END IF C C---- HISTORY "string" [OPTIONAL INPUT] C C History strings to be added to mtz o/p file HKLOUT C C ELSE IF (KEY.EQ.'HIST') THEN NHSOUT = NHSOUT + 1 IF (NHSOUT.GT.MAXHIS) THEN WRITE (LUNOUT,FMT=6020) MAXHIS 6020 FORMAT (' Sorry HISTORY now full no more lines accepted', + /,' ONLY ',I6,' Lines of HISTORY ALLOWED') NHSOUT = NHSOUT - 1 GO TO 10 END IF C C LXX = LEN(LINE) IF (LXX.GT.80) LXX = 80 HISOUT(NHSOUT) = LINE(IBEG(1) :LXX) KYHIST = .TRUE. GO TO 10 C C---- SYMMETRY C ======== C ELSE IF (KEY.EQ.'SYMM') THEN J = 2 C C ****************************************** CALL RDSYMM(J,LINE,IBEG,IEND,ITYP,FVALUE,NTOK,NAMSPG,NSPGRP, + PGNAME,NSYM,NSYMP,RSYM) C ******************************************** C KYSYMM = .TRUE. GO TO 10 C C C---- RESOlution C Resolution limits C ELSE IF (KEY.EQ.'RESO') THEN KYRESO = .TRUE. CALL RDRESO( $ 2,ITYP,FVALUE,NTOK,RMINWL,RMAXWL,SMINWL,SMAXWL) GO TO 10 C C---- RSCALE used for scaling. C Resolution limits - either 4sinsq/lamdasq or d C ELSE IF (KEY.EQ.'RSCA') THEN SCBLOG = .TRUE. CALL RDRESO( $ 2,ITYP,FVALUE,NTOK,RMINSC,RMAXSC,SMINSC,SMAXSC) GO TO 10 C C---- VPAT volume per atom - default = 10 C ELSE IF (KEY.EQ.'VPAT') THEN CALL GTPREA(2,VPA,NTOK,ITYP,FVALUE) GOTO 10 C C---- CONTents - unit cell content C ELSE IF (KEY.EQ.'CONT') THEN C C---- Check if Key_Word NRESidue given first, if so this C has priority and ignore keyword CONTents C IF (LOGRES) THEN WRITE (LUNOUT,6030) 6030 FORMAT (' You have already used the Key_Word NRESIDUE',/, + ' IGNORING this CONTENTS line') GO TO 10 END IF C C LOGCON = .TRUE. IAT = 0 C C DO 71 I=2,NTOK,2 IAT = IAT + 1 NW(IAT)= LINE( IBEG(I):(IBEG(I)+1) ) CALL CCPUPC(NW(IAT)) CALL GTPINT(I+1,NC(IAT),NTOK,ITYP,FVALUE) 71 CONTINUE C C NK=IAT GOTO 10 C C---- NRESIDUE Nres C C Nres iNrese number of residues expected in thein thein thein t C C A very approximate atom composition is calculated C C add adne ordered water per amino acid = cad = c C This is then taken as aken as aken as aken as aken as C as number of atoms in asymmetric unit. C ELSE IF (KEY.EQ.'NRES') THEN IF (NTOK.LE.1) THEN WRITE (LUNOUT,6040) 6040 FORMAT(' ERROR ---- for Key_Word NRESIDUES ',/, + ' NO value given for nres ... Stopping') CALL CCPERR(1,'Must give number of residues') END IF C C CALL GTPINT(2,NRES,NTOK,ITYP,FVALUE) NW(1) = 'C ' NC(1) = 5 * NRES NW(2) = 'N ' NC(2) = NINT(1.35 * NRES) NW(3) = 'O ' NC(3) = NINT(1.5 * NRES) NW(4) = 'H ' NC(4) = 8 * NRES NK = 4 LOGRES = .TRUE. GO TO 10 C C----BINS read number of ranges (bins), monitor interval C C Monitoring NUMBER _ program prints lots of information C about every nmon th reflection C C ELSE IF (KEY.EQ.'BINS') THEN C CALL GTPINT(2,NBIN,NTOK,ITYP,FVALUE) C CALL GTPINT(3,NMON,NTOK,ITYP,FVALUE) C GO TO 10 C C---- PLOT - turn ON of Fsq C ELSE IF (KEY.EQ.'PLOT') THEN IPLFSQ = 1 KEY = LINE(IBEG(2):IBEG(2)+3) CALL CCPUPC(KEY) IF(KEY.EQ.'OFF') THEN IPLFSQ = 0 ELSEIF(KEY.EQ.'ON') THEN IPLFSQ = 1 ENDIF GO TO 10 C c-harvest C C---- PRIVATE [priv] flag only C ======= C if present sets dir access for first C open of $HOME/DepositFiles to drwx------ C default is drwxr-xr-x C C ELSE IF (KEY .eq. 'PRIV') THEN PRIVATE = .true. GO TO 10 C C---- USECWD [usec] Flag only C ====== C if present file is opened in current working directory C default is o/p file to C $HOME/DepositFiles/ProjectName/DataSetName.ProgramName C ELSE IF (KEY .eq. 'USEC') THEN USECWD = .true. GO TO 10 C C---- PROJECTNAME [pname] C =========== C C if given with DATASET then harvest will o/p a file C this project_name should be always used for the C one structure determination C no default C ELSE IF (KEY .eq. 'PNAM') THEN GOT_PNAME = .true. PNAME_KEYWRD = LINE(IBEG(2) :IEND(2)) GO TO 10 C C---- DATASETNAME [dname] C =========== C C if given with PROJECT then harvest will o/p a file C this dataset_name is the name of one of the diffraction C data sets used in a particular project C no default C ELSE IF (KEY .eq. 'DNAM') THEN DNAME_KEYWRD = LINE(IBEG(2) :IEND(2)) GOT_DNAME = .true. GO TO 10 C C---- ROWLIMIT [rsize] C ======== C C if given sets the number of characters in a row C for the output file, default is 132 C TRUE mmcif default is 80 C ELSE IF (KEY .eq. 'RSIZ') THEN NROWLIMIT = 80 IF (NTOK.GE.2) CALL GTPINT(2,NROWLIMIT,NTOK,ITYPE,FVALUE) GO TO 10 C C---- NOHARVEST Flag only C ====== C disable writing of harvest files C ELSE IF (KEY .eq. 'NOHA') THEN NOHARVEST = .true. GO TO 10 c-harvest C---- SCALE WILSON | C set scale factor C WILSON [default] use scale from Wilson plot C ELSE IF (KEY .EQ. 'SCAL') THEN IF (NTOK .GT. 1) THEN IF (ITYP(2) .EQ. 1) THEN CALL CCPUPC(CVALUE(2)) IF (CVALUE(2) .EQ. 'WILS') THEN LAPSCL = .TRUE. ENDIF ELSEIF(ITYP(2) .EQ. 2) THEN CALL GTPREA(2,SCALEF,NTOK,ITYP,FVALUE) LAPSCL = .FALSE. ENDIF ENDIF GO TO 10 C ELSEIF (KEY.EQ.'END') THEN GOTO 50 ELSE WRITE (LUNOUT,FMT=6022) LINE(1:LENSTR(LINE)) 6022 FORMAT (' Error: Key_Word line NOT Understood:-',/,' ',A) GO TO 10 END IF C C 50 CONTINUE C C---- Check if WILSON Key_Words CONT or NRES used C IF (.NOT.LOGCON .AND. .NOT.LOGRES) THEN WRITE (LUNOUT,6026) 6026 FORMAT ( $ ' No CONTENTS or NRESIDUE Key_Word given with WILSON',/, + ' Number of residues will be estimated from volume', + ' assuming 50% solvent.') ENDIF C C END C C C SUBROUTINE LABEL(VMIN,VMAX,LOWV,INCREM,NTICK,IPOWER) C ==================================================== 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 C POWER = LOG10(VMAX-VMIN) C C---- Make sure to correct for fact that int 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 C C SUBROUTINE LSQ(X,Y,NOBS,IOBS,SIG,MWT,A,B,SIGA,SIGB,CHI2,Q) C ===================================================== C C 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 very general lsq program C From "Numerical recipes" C C date January 1988 C Given a set of points x(i), y(i) with SD sig(i), i=iobs to nobs, C fit 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 implicit none C C INTEGER I,MWT,NOBS REAL A,B,SIGA,SIGB,CHI2,Q,SS,WT,SIGDAT REAL X(NOBS),Y(NOBS),SIG(NOBS),T REAL SXX,SYY,SXY,SX,SY,ST2,SXOSS,TMP LOGICAL EMPTY_SHELL C SXX=0 SYY=0 SXY=0 SX=0 SY=0 ST2=0 B=0 EMPTY_SHELL=.FALSE. C IF (MWT.EQ.0)THEN DO 10 I=1,NOBS SIG(I)=1.0 10 CONTINUE ENDIF C C SS=0 DO 11 I=IOBS,NOBS IF(SIG(I).EQ.0.0) GO TO 11 WT=1.0/(SIG(I)**2) SS=SS+WT SXX=SXX+ X(I)*X(I)* WT SX=SX+X(I)*WT SXY=SXY+X(I)*Y(I)*WT SYY=SYY+Y(I)*Y(I)*WT SY=SY+Y(I)*WT 11 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 x*x*wt C ---- syy= sum of y*Y*wt C ---- sxy= sum of x*y*wt C C Find a and b following recipe p 508 C SXOSS=SX/SS C C DO 13 I=IOBS,NOBS IF(SIG(I).EQ.0.0 ) THEN EMPTY_SHELL=.TRUE. GO TO 13 ENDIF T=(X(I)-SXOSS)/SIG(I) ST2=ST2+T*T B=B+T*Y(I)/SIG(I) 13 CONTINUE C C B=B/ST2 A=(SY-SX*B)/SS SIGA=SQRT((1.0+SX*SX/(SS*ST2))/SS) SIGB=SQRT(1.0/ST2) C C --- Q is returned as the correlation coefficient C C This is TJO original...... 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 16 I=IOBS,NOBS IF(SIG(I).EQ.0.0 ) GO TO 16 CHI2=CHI2+((Y(I)-A-B*X(I))/SIG(I))**2 TMP=(Y(I)-A-B*X(I))/SIG(I) SIG(I)=ABS(TMP) 16 CONTINUE C 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 ************* the above is not implemented ***** C IF (EMPTY_SHELL) THEN CALL CCPERR(3, 'WARNING: Check empty shell of data with hklview') ENDIF RETURN END C C SUBROUTINE NEGIAS(XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ) C ================================================== C C IMPLICIT NONE C C C C---- Subroutine truncates a normal distribution N(XI-SDI/SIGMA, C SDI zero and returns the mean and s.d. of the positive part C together with C the mean and s.d. of the square root of the positive part. C C---- Essentially this routine gives the means and s.d's of the post C distributions for the true intensity and structure factor modu C acentric reflections where an intensity of xi has been observe C with s.d. sdi and the wilsons distribution for the region of C reciprocal space in which the reflection lies has variance sig C C---- Accuracy - better than 5 percent (i think) C C---- This version uses cubic spline interpolation C C---- Arguments: C XI - Mean of distribution to be truncated. C SDI - S.D. of distribution to be truncated. C SIGMA - Variance of Wilson's distribution C XJ - Mean of truncated positive part. C SDJ - S.D. of truncated positive part. C XF - Mean of square root of positive part. C SDF - S.D. of square root of positive part. C IREJ - ON ENTRY: = -1 if only XF and SDF required. C = 0 if all XJ,SDJ,XF, and SDJ required. C = 1 if only XJ and SDJ required. C ON EXIT: = -1 if (XI-SDI/SIGMA) < -4.0*SDI C = 0 if all O.K. C = 1 if observation is essentially imposs C (i.e. XI < -3.7*SDI) C C C Arrays ZJ,ZJSD,ZF,ZFSD contain the tabulated values of XJ,SDJ, C for XI = -4.0(0.1)3.0. C Linear interpolation is used in these tables. C C C .. Scalar Arguments .. REAL SDF,SDI,SDJ,SIGMA,XF,XI,XJ INTEGER IREJ C .. C .. Local Scalars .. REAL H INTEGER IOP C .. C .. Local Arrays .. REAL ZF(71),ZFSD(71),ZJ(71),ZJSD(71) REAL ZFD2(71),ZFSDD2(71),ZJD2(71),ZJSDD2(71),RANGES(71) REAL ZFIP,ZFSDIP,ZJIP,ZJSDIP C .. C .. External Subroutines EXTERNAL SPLINT C .. C .. Intrinsic Functions .. INTRINSIC SQRT C .. C .. Data statements .. DATA ZJ/0.226,0.230,0.235,0.240,0.246,0.251,0.257,0.263,0.270, + 0.276,0.283,0.290,0.298,0.306,0.314,0.323,0.332,0.341,0.351, + 0.362,0.373,0.385,0.397,0.410,0.424,0.439,0.454,0.470,0.487, + 0.505,0.525,0.545,0.567,0.590,0.615,0.641,0.668,0.698,0.729, + 0.762,0.798,0.835,0.875,0.917,0.962,1.009,1.059,1.112,1.167, + 1.226,1.287,1.352,1.419,1.490,1.563,1.639,1.717,1.798,1.882, + 1.967,2.055,2.145,2.236,2.329,2.422,2.518,2.614,2.710,2.808, + 2.906,3.004/ DATA ZJSD/0.217,0.221,0.226,0.230,0.235,0.240,0.245,0.250,0.255, + 0.261,0.267,0.273,0.279,0.286,0.292,0.299,0.307,0.314,0.322, + 0.330,0.339,0.348,0.357,0.367,0.377,0.387,0.398,0.409,0.421, + 0.433,0.446,0.459,0.473,0.488,0.503,0.518,0.535,0.551,0.568, + 0.586,0.604,0.622,0.641,0.660,0.679,0.698,0.718,0.737,0.757, + 0.776,0.795,0.813,0.831,0.848,0.865,0.881,0.895,0.909,0.921, + 0.933,0.943,0.953,0.961,0.968,0.974,0.980,0.984,0.988,0.991, + 0.994,0.996/ DATA ZF/0.423,0.428,0.432,0.437,0.442,0.447,0.453,0.458,0.464, + 0.469,0.475,0.482,0.488,0.495,0.502,0.509,0.516,0.524,0.532, + 0.540,0.549,0.557,0.567,0.576,0.586,0.597,0.608,0.619,0.631, + 0.643,0.656,0.670,0.684,0.699,0.714,0.730,0.747,0.765,0.783, + 0.802,0.822,0.843,0.865,0.887,0.911,0.935,0.960,0.987,1.014, + 1.042,1.070,1.100,1.130,1.161,1.192,1.224,1.257,1.289,1.322, + 1.355,1.388,1.421,1.454,1.487,1.519,1.551,1.583,1.615,1.646, + 1.676,1.706/ DATA ZFSD/0.216,0.218,0.220,0.222,0.224,0.226,0.229,0.231,0.234, + 0.236,0.239,0.241,0.244,0.247,0.250,0.253,0.256,0.259,0.262, + 0.266,0.269,0.272,0.276,0.279,0.283,0.287,0.291,0.295,0.298, + 0.302,0.307,0.311,0.315,0.319,0.324,0.328,0.332,0.337,0.341, + 0.345,0.349,0.353,0.357,0.360,0.364,0.367,0.369,0.372,0.374, + 0.375,0.376,0.377,0.377,0.377,0.376,0.374,0.372,0.369,0.366, + 0.362,0.358,0.353,0.348,0.343,0.338,0.332,0.327,0.321,0.315, + 0.310,0.304/ DATA RANGES/-4.0,-3.9,-3.8,-3.7,-3.6,-3.5,-3.4,-3.3,-3.2,-3.1, + -3.0,-2.9,-2.8,-2.7,-2.6,-2.5,-2.4,-2.3,-2.2,-2.1, + -2.0,-1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1, + -1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1, + 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, + 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, + 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, + 3.0/ C .. C Second derivatives DATA ZJD2/ 0.000000, 0.001763,-0.001053, 0.002447,-0.002736, + 0.002496,-0.001248, 0.002496,-0.002738, 0.002455,-0.001080, + 0.001867,-0.000386,-0.000321, 0.001672,-0.000366,-0.000209, + 0.001200, 0.001408,-0.000833, 0.001924,-0.000863, 0.001530, + 0.000743, 0.001497,-0.000730, 0.001422, 0.001042, 0.000409, + 0.003320,-0.001690, 0.003439,-0.000067, 0.002831, 0.000744, + 0.000194, 0.004481,-0.000117, 0.001988, 0.004166,-0.000651, + 0.004439, 0.000894, 0.003984, 0.001171, 0.003333, 0.003499, + 0.000671, 0.005819, 0.000051, 0.005977, 0.000040, 0.005861, + 0.000515, 0.004080, 0.001166, 0.003257, 0.003807,-0.000487, + 0.004139, 0.001931, 0.000138, 0.003520,-0.002219, 0.005356, + -0.001204,-0.000539, 0.003359,-0.000896, 0.000224, 0.000000/ DATA ZJSDD2/ 0.000000, 0.002154,-0.002615, 0.002306,-0.000610, + 0.000132, 0.000081,-0.000455, 0.001741,-0.000507, 0.000288, + -0.000643, 0.002285,-0.002495, 0.001697, 0.001707,-0.002526, + 0.002398,-0.001065, 0.001863,-0.000386,-0.000318, 0.001660, + -0.000320,-0.000378, 0.001833,-0.000955, 0.001985,-0.000986, + 0.001958,-0.000845, 0.001423, 0.001155,-0.000042,-0.000989, + 0.003996,-0.002996, 0.001988, 0.001044,-0.000164,-0.000390, + 0.001723,-0.000504, 0.000292,-0.000664, 0.002365,-0.002797, + 0.002821,-0.002488, 0.001129,-0.002029, 0.000987,-0.001921, + 0.000697,-0.000866,-0.003232, 0.001793,-0.003942, 0.001974, + -0.003954, 0.001841,-0.003410,-0.000199,-0.001794, 0.001375, + -0.003707, 0.001452,-0.002103, 0.000961,-0.001740, 0.000000/ DATA ZFD2/ 0.000000,-0.002027, 0.002109,-0.000408,-0.000478, + 0.002318,-0.002795, 0.002861,-0.002650, 0.001740, 0.001688, + -0.002493, 0.002284,-0.000642, 0.000284,-0.000495, 0.001696, + -0.000289,-0.000541, 0.002451,-0.003264, 0.004606,-0.003160, + 0.002032, 0.001031,-0.000155,-0.000411, 0.001800,-0.000789, + 0.001356, 0.001367,-0.000822, 0.001921,-0.000863, 0.001530, + 0.000743, 0.001499,-0.000738, 0.001455, 0.000920, 0.000868, + 0.001610,-0.001307, 0.003616,-0.001159, 0.001021, 0.003074, + -0.001319, 0.002200,-0.001484, 0.003734,-0.001451, 0.002071, + -0.000831, 0.001253, 0.001818,-0.002525, 0.002284,-0.000610, + 0.000155,-0.000011,-0.000113, 0.000463,-0.001738, 0.000488, + -0.000215, 0.000374,-0.001279,-0.001259, 0.000315, 0.000000/ DATA ZFSDD2/ 0.000000, 0.000011,-0.000045, 0.000170,-0.000634, + 0.002366,-0.002829, 0.002951,-0.002976, 0.002951,-0.002829, + 0.002366,-0.000633, 0.000167,-0.000034,-0.000029, 0.000151, + -0.000576, 0.002151,-0.002030,-0.000033, 0.002162,-0.002615, + 0.002299,-0.000579, 0.000019, 0.000501,-0.002026, 0.001601, + 0.001620,-0.002080, 0.000701,-0.000724, 0.002193,-0.002049, + 0.000003, 0.002037,-0.002151, 0.000565,-0.000111,-0.000122, + 0.000599,-0.002273, 0.002495,-0.001707,-0.001666, 0.002372, + -0.001821,-0.001087, 0.000170, 0.000407,-0.001799, 0.000788, + -0.001355,-0.001369, 0.000831,-0.001954, 0.000986,-0.001988, + 0.000967,-0.001878, 0.000544,-0.000300, 0.000655,-0.002319, + 0.002621,-0.002166, 0.000045, 0.001988,-0.001997, 0.000000/ C C IOP = IREJ IREJ = 0 C C---- Is observation possible? C H = XI/SDI IF (H+3.7.GT.0) THEN C C---- Locate range of distribution C H = H - SDI/SIGMA IF (H+4.0.LE.0) THEN C C---- Distribution is essentially negative. C IREJ = -1 C ELSE IF (H-3.0.LT.0) THEN C C---- Interpolate in the tables C IF (IOP.GE.0) THEN CALL SPLINT(RANGES,ZJ,ZJD2,71,H,ZJIP) XJ = ZJIP*SDI CALL SPLINT(RANGES,ZJSD,ZJSDD2,71,H,ZJSDIP) SDJ = ZJSDIP*SDI END IF C C IF (IOP.LE.0) THEN SDF = SQRT(SDI) CALL SPLINT(RANGES,ZF,ZFD2,71,H,ZFIP) XF = ZFIP*SDF CALL SPLINT(RANGES,ZFSD,ZFSDD2,71,H,ZFSDIP) SDF = ZFSDIP*SDF END IF C C C RETURN ELSE C C---- Distribution is essentially all positive C H = H*SDI C C IF (IOP.GE.0) THEN XJ = H SDJ = SDI END IF C C IF (IOP.LE.0) THEN XF = SQRT(H) SDF = SDI*0.5/XF END IF C C RETURN END IF C C ELSE IREJ = 1 END IF C C XJ = 0.0 SDJ = 1.0 XF = 0.0 SDF = 1.0 C C END C C C SUBROUTINE NEGICS(XI,SDI,SIGMA,XJ,SDJ,XF,SDF,IREJ) C ================================================== C C IMPLICIT NONE C C C C---- Subroutine truncates a normal distribution N(XI-SDI/SIGMA, C SDI zero and returns the mean and s.d. of the positive part C together with C the mean and s.d. of the square root of the positive part. C C---- Essentially this routine gives the means and s.d's of the post C distributions for the true intensity and structure factor modu C centric reflections where an intensity of xi has been observed C with s.d. sdi and the wilsons distribution for the region of C reciprocal space in which the reflection lies has variance sig C C---- Accuracy - better than 5 percent (i think) C C---- This version uses cubic spline interpolation C C---- Arguments: C XI - Mean of distribution to be truncated. C SDI - S.D. of distribution to be truncated. C SIGMA - Variance of Wilson's distribution C XJ - Mean of truncated positive part. C SDJ - S.D. of truncated positive part. C XF - Mean of square root of positive part. C SDF - S.D. of square root of positive part. C IREJ - ON ENTRY: = -1 if only XF and SDF required. C = 0 if all XJ,SDJ,XF, and SDJ required. C = 1 if only XJ and SDJ required. C ON EXIT: = -1 if (XI-SDI/SIGMA) < -4.0*SDI C = 0 if all O.K. C = 1 if observation is essentially imposs C (i.e. XI < -3.7*SDI) C C C Arrays ZJ,ZJSD,ZF,ZFSD contain the tabulated values of XJ,SDJ, C for XI = -4.0(0.1)4.0. C linear interpolation is used in these tables. C C C .. Scalar Arguments .. REAL SDF,SDI,SDJ,SIGMA,XF,XI,XJ INTEGER IREJ C .. C .. Local Scalars .. REAL H,XI2,XI4,XSDF,XXF INTEGER IOP C .. C .. Local Arrays .. REAL ZF(81),ZFSD(81),ZJ(81),ZJSD(81) REAL ZFD2(81),ZFSDD2(81),ZJD2(81),ZJSDD2(81),RANGES(81) REAL ZFIP,ZFSDIP,ZJIP,ZJSDIP C .. C .. External Subroutines EXTERNAL SPLINT C .. C .. Intrinsic Functions .. INTRINSIC SQRT C .. C .. Data statements .. DATA ZJ/0.114,0.116,0.119,0.122,0.124,0.127,0.130,0.134,0.137, + 0.141,0.145,0.148,0.153,0.157,0.162,0.166,0.172,0.177,0.183, + 0.189,0.195,0.202,0.209,0.217,0.225,0.234,0.243,0.253,0.263, + 0.275,0.287,0.300,0.314,0.329,0.345,0.363,0.382,0.402,0.425, + 0.449,0.475,0.503,0.534,0.567,0.603,0.642,0.684,0.730,0.779, + 0.833,0.890,0.952,1.018,1.089,1.164,1.244,1.327,1.416,1.508, + 1.603,1.703,1.805,1.909,2.015,2.123,2.233,2.343,2.453,2.564, + 2.674,2.784,2.894,3.003,3.112,3.220,3.328,3.435,3.541,3.647, + 3.753,3.962/ DATA ZJSD/0.158,0.161,0.165,0.168,0.172,0.176,0.179,0.184,0.188, + 0.192,0.197,0.202,0.207,0.212,0.218,0.224,0.230,0.236,0.243, + 0.250,0.257,0.265,0.273,0.282,0.291,0.300,0.310,0.321,0.332, + 0.343,0.355,0.368,0.382,0.397,0.412,0.428,0.445,0.463,0.481, + 0.501,0.521,0.543,0.565,0.589,0.613,0.638,0.664,0.691,0.718, + 0.745,0.773,0.801,0.828,0.855,0.881,0.906,0.929,0.951,0.971, + 0.989,1.004,1.018,1.029,1.038,1.044,1.049,1.052,1.054,1.054, + 1.053,1.051,1.049,1.047,1.044,1.041,1.039,1.036,1.034,1.031, + 1.029,1.028/ DATA ZF/0.269,0.272,0.276,0.279,0.282,0.286,0.289,0.293,0.297, + 0.301,0.305,0.309,0.314,0.318,0.323,0.328,0.333,0.339,0.344, + 0.350,0.356,0.363,0.370,0.377,0.384,0.392,0.400,0.409,0.418, + 0.427,0.438,0.448,0.460,0.471,0.484,0.498,0.512,0.527,0.543, + 0.560,0.578,0.597,0.618,0.639,0.662,0.687,0.713,0.740,0.769, + 0.800,0.832,0.866,0.901,0.938,0.976,1.016,1.057,1.098,1.140, + 1.183,1.227,1.270,1.313,1.356,1.398,1.439,1.480,1.519,1.558, + 1.595,1.632,1.667,1.701,1.735,1.767,1.799,1.829,1.859,1.889, + 1.917,1.945/ DATA ZFSD/0.203,0.205,0.207,0.209,0.211,0.214,0.216,0.219,0.222, + 0.224,0.227,0.230,0.233,0.236,0.239,0.243,0.246,0.250,0.253, + 0.257,0.261,0.265,0.269,0.273,0.278,0.283,0.288,0.293,0.298, + 0.303,0.309,0.314,0.320,0.327,0.333,0.340,0.346,0.353,0.361, + 0.368,0.375,0.383,0.390,0.398,0.405,0.413,0.420,0.427,0.433, + 0.440,0.445,0.450,0.454,0.457,0.459,0.460,0.460,0.458,0.455, + 0.451,0.445,0.438,0.431,0.422,0.412,0.402,0.392,0.381,0.370, + 0.360,0.349,0.339,0.330,0.321,0.312,0.304,0.297,0.290,0.284, + 0.278,0.272/ DATA RANGES/-4.0,-3.9,-3.8,-3.7,-3.6,-3.5,-3.4,-3.3,-3.2,-3.1, + -3.0,-2.9,-2.8,-2.7,-2.6,-2.5,-2.4,-2.3,-2.2,-2.1, + -2.0,-1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1, + -1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1, + 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, + 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, + 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, + 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, + 4.0/ C C..Second derivatives DATA ZJD2/ 0.000000, 0.001458, 0.000166,-0.002124, 0.002331, + -0.001198, 0.002461,-0.002647, 0.002128, 0.000134,-0.002665, + 0.004526,-0.003439, 0.003230,-0.003480, 0.004692,-0.003287, + 0.002455,-0.000533,-0.000321, 0.001818,-0.000951, 0.001988, + -0.000999, 0.002008,-0.001032, 0.002122,-0.001455, 0.003697, + -0.001332, 0.001633, 0.000801, 0.001164, 0.000542, 0.002667, + 0.000791, 0.000168, 0.004538,-0.000320, 0.002744, 0.001346, + 0.003871, 0.001168, 0.003455, 0.003011, 0.002500, 0.004988, + 0.001548, 0.006820, 0.001172, 0.006493, 0.002857, 0.006078, + 0.002830, 0.006602, 0.000761, 0.008354, 0.001820, 0.002364, + 0.006724, 0.000737, 0.002326, 0.001959, 0.001839, 0.002684, + -0.000574,-0.000389, 0.002133,-0.002145, 0.000451, 0.000342, + -0.001817, 0.000925,-0.001883, 0.000608,-0.000548,-0.004418, + 0.012220,-0.044459, 0.165615, 0.000000/ DATA ZJSDD2/ 0.000000, 0.002141,-0.002563, 0.002109, 0.000125, + -0.002610, 0.004314,-0.002648, 0.000276, 0.001542,-0.000444, + 0.000235,-0.000496, 0.001749,-0.000499, 0.000247,-0.000488, + 0.001705,-0.000333,-0.000372, 0.001821,-0.000913, 0.001829, + -0.000405,-0.000209, 0.001239, 0.001251,-0.000245,-0.000272, + 0.001332, 0.000942, 0.000899, 0.001461,-0.000743, 0.001509, + 0.000706, 0.001667,-0.001375, 0.003832,-0.001952, 0.003976, + -0.001952, 0.003831,-0.001372, 0.001658, 0.000740, 0.001381, + -0.000263,-0.000327, 0.001573, 0.000036,-0.001717, 0.000833, + -0.001614,-0.000377,-0.002877,-0.000116,-0.002660,-0.001243, + -0.004368, 0.000713,-0.004483,-0.000780,-0.004398, 0.000373, + -0.003092,-0.000003,-0.002897,-0.000408,-0.001470, 0.000290, + 0.000311,-0.001534,-0.000176, 0.002238,-0.002777, 0.002869, + -0.002697, 0.001919, 0.001020, 0.000000/ DATA ZFD2/ 0.000000, 0.001997,-0.001988,-0.000044, 0.002166, + -0.002618, 0.002306,-0.000607, 0.000121, 0.000122,-0.000609, + 0.002314,-0.002648, 0.002276,-0.000455,-0.000458, 0.002285, + -0.002682, 0.002442,-0.001087, 0.001908,-0.000544, 0.000269, + -0.000532, 0.001860,-0.000909, 0.001776,-0.000194,-0.000999, + 0.004190,-0.003760, 0.004852,-0.003647, 0.003738, 0.000697, + -0.000526, 0.001406, 0.000902, 0.000985, 0.001156, 0.000389, + 0.003287,-0.001537, 0.002863, 0.002087, 0.000791, 0.000750, + 0.002211, 0.002408, 0.000159, 0.002956, 0.000016, 0.002978, + 0.000070, 0.002740, 0.000968,-0.000614, 0.001487, 0.000665, + 0.001855,-0.002083, 0.000479, 0.000169,-0.001156,-0.001547, + 0.001342,-0.003820, 0.001937,-0.003929, 0.001781,-0.003195, + -0.000999, 0.001190,-0.003761, 0.001853,-0.003650, 0.000748, + 0.000657,-0.003375, 0.000844, 0.000000/ DATA ZFSDD2/ 0.000000,-0.000041, 0.000165,-0.000619, 0.002309, + -0.002618, 0.002163,-0.000034,-0.002029, 0.002148,-0.000564, + 0.000109, 0.000128,-0.000622, 0.002360,-0.002817, 0.002909, + -0.002817, 0.002360,-0.000624, 0.000136, 0.000079,-0.000452, + 0.001730,-0.000466, 0.000136,-0.000077, 0.000172,-0.000612, + 0.002276,-0.002491, 0.001688, 0.001738,-0.002640, 0.002823, + -0.002652, 0.001783, 0.001518,-0.001857,-0.000091, 0.002220, + -0.002789, 0.002935,-0.002952, 0.002873,-0.002541, 0.001291, + -0.002622, 0.003196,-0.004163, 0.001455,-0.001655,-0.000833, + -0.001014,-0.001110,-0.000547,-0.002701,-0.000649,-0.000702, + -0.002543,-0.001127, 0.001050,-0.003073,-0.000756, 0.000098, + 0.000364,-0.001553,-0.000153, 0.002164,-0.002505, 0.001854, + 0.001089,-0.000210,-0.000251, 0.001213, 0.001398,-0.000805, + 0.001823,-0.000486, 0.000122, 0.000000/ C C IOP = IREJ IREJ = 0 C C---- Is observation possible? C H = XI/SDI C C IF (H+3.7.GT.0) THEN H = H - SDI/ (SIGMA*2.0) C C IF (H+4.0.LE.0) THEN IREJ = -1 ELSE IF (H-4.0.LT.0) THEN C C---- Interpolate in the tables C IF (IOP.GE.0) THEN CALL SPLINT(RANGES,ZJ,ZJD2,71,H,ZJIP) XJ = ZJIP*SDI CALL SPLINT(RANGES,ZJSD,ZJSDD2,71,H,ZJSDIP) SDJ = ZJSDIP*SDI END IF C C IF (IOP.LE.0) THEN SDF = SQRT(SDI) CALL SPLINT(RANGES,ZF,ZFD2,71,H,ZFIP) XF = ZFIP*SDF CALL SPLINT(RANGES,ZFSD,ZFSDD2,71,H,ZFSDIP) SDF = ZFSDIP*SDF END IF C C RETURN ELSE C C---- Distribution is essentially all positive C XI2 = 1.0/ (H*H) XI4 = XI2*XI2 XXF = (1.0-0.375*XI2-87.0/128.0*XI4)*SQRT(H) XSDF = SQRT((30.0/64.0*XI4+0.25*XI2)*H) C C IF (IOP.GE.0) THEN XJ = H*SDI* (1.0-0.5*XI2-0.75*XI4) SDJ = SDI*XSDF*2.0*XXF END IF C C IF (IOP.LE.0) THEN SDF = SQRT(SDI) XF = XXF*SDF SDF = XSDF*SDF END IF C C RETURN END IF C C ELSE IREJ = 1 END IF C C XJ = 0.0 SDJ = 1.0 XF = 0.0 SDF = 1.0 C C END C C C SUBROUTINE PLOTR(XY,NXY,ILSQ,CXY,TITLEX,TITLEY,LEGEND) C ================================================= 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 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,ILSQ CHARACTER CXY*1 C .. C .. Array Arguments .. REAL XY(2,100) CHARACTER LEGEND*60,TITLEX*20,TITLEY*20 C .. C .. Local Scalars .. REAL XMAX,XMIN,XRANGE,YMAX,YMIN, + YRANGE INTEGER I,I1,I2,INCX,INCY,IPOWRX,IPOWRY, + ISCALE,J,K,KINT,LOWX,LOWY,NOBS, + NTICKX,NTICKY CHARACTER BAR*1,BLANK*1,CORNER*1,DASH*1,ONE*1,STAR*1, + TICK*1 C .. C .. Local Arrays .. CHARACTER KCHAR(5)*1,PLOT(129,60)*1 C .. C .. External Subroutines .. EXTERNAL ITRANS,LABEL C .. C .. Intrinsic Functions .. INTRINSIC INT,SQRT C .. C .. Data statements .. C C---- Choose most suitable characters available on printer. C DATA BLANK/' '/,BAR/'|'/,DASH/'-'/,ONE/'1'/, + CORNER/'+'/ DATA TICK/'+'/,STAR/'*'/ C .. C C Try 0 0 - not set ???????????????????????????????? LOWX = 0 LOWY = 0 C C NOBS = NXY IF (NOBS.GT.60) THEN WRITE(6,'(a,i3,a)') 1 ' ** RESET NOBS TO 60 FOR PLOTTING ** ',NXY,' OBSERVATIONS **' NOBS = 60 END IF C C XMIN = 1.0E30 XMAX = -1.0E30 YMIN = 1.0E30 YMAX = -1.0E30 C C---- Determine outer bounds of x and y values C IF (NOBS.GE.1) THEN DO 50 I = 1,NOBS 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 C CALL LABEL(XMIN,XMAX,LOWX,INCX,NTICKX,IPOWRX) CALL LABEL(YMIN,YMAX,LOWY,INCY,NTICKY,IPOWRY) 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 C C C PLOT(1,1) = ONE PLOT(9,55) = CORNER DO 80 I = 1,60 I1 = I + 40 I2 = I + 20 PLOT(I1,58) = ' ' IF(I2.LE.60)PLOT(2,I2) = ' ' PLOT(I+40,60) = LEGEND(I:I) C IF(I.GT.20) GO TO 80 PLOT(I1,58) = TITLEX(I:I) PLOT(2,I2) = TITLEY(I:I) 80 CONTINUE C C---- If x axis labels can fit into 5 characters, label x axis. C 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 CALL ITRANS(KINT,KCHAR) DO 90 I2 = 1,5 IF (KCHAR(I2).NE.BLANK) THEN PLOT(J,56) = KCHAR(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 CALL ITRANS(ISCALE,KCHAR) PLOT(62,58) = STAR PLOT(63,58) = ONE PLOT(64,58) = '0' PLOT(65,58) = '*' PLOT(66,58) = '*' I = 67 DO 110 J = 1,5 IF (KCHAR(J).NE.BLANK) THEN PLOT(I,58) = KCHAR(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 CALL ITRANS(KINT,KCHAR) DO 120 I2 = 1,5 K = I2 + 3 PLOT(K,J) = KCHAR(I2) 120 CONTINUE END IF 130 CONTINUE C C---- Indicate multiple of 10. C IF (IPOWRY.NE.0) THEN ISCALE = -IPOWRY CALL ITRANS(ISCALE,KCHAR) PLOT(2,42) = STAR PLOT(2,43) = ONE PLOT(2,44) = '0' PLOT(2,45) = '*' PLOT(2,46) = '*' I = 47 DO 140 J = 1,5 IF (KCHAR(J).NE.BLANK) THEN PLOT(2,I) = KCHAR(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) IF(J.GE.1 .AND. J.LE.129 .AND. 1 K.GE.1 .AND. K.LE.60) 1 PLOT(J,K) = CXY 150 CONTINUE C C---- Plot is ready for printer. C WRITE (6,FMT=9020) PLOT END IF END IF RETURN C C---- Format statements C 9000 FORMAT (' EXCLUDE',I4,3F10.4) 9010 FORMAT (I12,3F10.4) 9020 FORMAT (59 (129A1,/),129A1) END C C C SUBROUTINE SETCLL(C) C ==================== C C IMPLICIT NONE C C C C---- Calculate reciprocal cell dimensions R from C direct cell dimensions C C C C C C C .. Array Arguments .. REAL C(6) C .. C .. Arrays in Common .. REAL R C .. C .. Local Scalars .. REAL SUM,V INTEGER JDO10,JDO20,JDO30 C .. C .. Intrinsic Functions .. INTRINSIC ABS,COS,SIN,SQRT C .. C .. Common blocks .. COMMON /CELLDM/R(6) C .. SAVE /CELLDM/ C C SUM = 0.0 C C DO 10 JDO10 = 4,6 IF (ABS(C(JDO10)).LT.0.1) C(JDO10) = 90.0 C(JDO10) = C(JDO10)*0.0174532 SUM = C(JDO10)*0.5 + SUM 10 CONTINUE C C V = 1.0 C C DO 20 JDO20 = 4,6 V = SIN(SUM-C(JDO20))*V 20 CONTINUE C C V = C(1)*2.0*C(2)*C(3)*SQRT(SIN(SUM)*V) R(1) = C(2)*C(3)*SIN(C(4))/V R(2) = C(1)*C(3)*SIN(C(5))/V R(3) = C(1)*C(2)*SIN(C(6))/V R(4) = (COS(C(5))*COS(C(6))-COS(C(4)))/ (SIN(C(5))*SIN(C(6))) R(5) = (COS(C(6))*COS(C(4))-COS(C(5)))/ (SIN(C(6))*SIN(C(4))) R(6) = (COS(C(4))*COS(C(5))-COS(C(6)))/ (SIN(C(4))*SIN(C(5))) C C DO 30 JDO30 = 4,6 C(JDO30) = C(JDO30)/0.0174532 30 CONTINUE C C R(4) = R(4)*2.0*R(2)*R(3) R(5) = R(5)*2.0*R(1)*R(3) R(6) = R(6)*2.0*R(1)*R(2) C C END C C C REAL FUNCTION STHLQ4(IN) C ======================== C C IMPLICIT NONE C C C---- Calculate 4(sintheta/lambda)**2 C C C C .. Array Arguments .. INTEGER IN(3) C .. C .. Arrays in Common .. REAL R C .. C .. Local Scalars .. INTEGER JDO10 C .. C .. Local Arrays .. REAL A(3) C .. C .. Common blocks .. COMMON /CELLDM/R(6) C .. SAVE /CELLDM/ C C STHLQ4 = 0.0 C C DO 10 JDO10 = 1,3 A(JDO10) = IN(JDO10) STHLQ4 = (A(JDO10)*R(JDO10))**2 + STHLQ4 10 CONTINUE C C STHLQ4 = R(4)*A(2)*A(3) + STHLQ4 + R(5)*A(1)*A(3) + R(6)*A(1)*A(2) C C END C C C SUBROUTINE WILPLT(NSYM,NSYMP,NREFWL,NRANGE,WILSC,WILBT,IWITCH) C ====================================================== C C**********************ERRORvariable CORRN Not defined C setting this to 1.0 as I dont know C what it means C INTEGER IWITCH,NSYM,NSYMP,NREFWL,IFFAV1,IFFAV2,NOPT,I,NRFF,JJ REAL WILSC(2),WILBT(2),FF1(100),FF2(100),XLSQ(100),YLSQ(100), + SIG(100),FFAVV1(100),FFAVV2(100),FFAV1,FFAV2,RESI, + SIGRMS,A,B,CHI2,Q,SIGA,SIGB,CORRN CHARACTER LEGEND*60,TITLEX*20,TITLEY*20 C INTEGER NBIN,NMON,IPLFSQ,NK REAL RMAX,RMIN,SMIN,SMAX,SMINSC,SMAXSC,RMINSC,RMAXSC, $ VPA REAL SI,SN,SR,SW,XY,XYAV,FFS,PI,DELDSS,F000, + VOL,RR,CELLAS,RVOL,FRAC INTEGER NRFUNQ,NATMS,NELEC,ISMAX C C COMMON /WILKEY/ NBIN,NMON,RMAX,RMIN,SMIN,SMAX,SMINSC, + SMAXSC,RMINSC,RMAXSC,VPA, + IPLFSQ,NK COMMON /PREWIL/ SI(100),SN(100),SR(100),SW(100), + NRFUNQ(100),XY(2,100),XYAV(2,100), + FFS(500),PI,DELDSS,NATMS,NELEC,F000, + ISMAX,VOL,RR(3,3,6),CELLAS(6),RVOL,FRAC CHARACTER NW*4 COMMON /WILCHR/ NW(20) INTEGER NC,KatmWghts REAL FormFactAcoeff,FormFactBcoeff COMMON /WILFRM/ NC(20),KatmWghts(20),FormFactAcoeff(5,20), + FormFactBcoeff(4,20) DATA SIG/100*0.0/ DATA XLSQ/100*0.0/ DATA YLSQ/100*0.0/ C C---- ????????????? CORRN = 1.0 C---- ????????????? C C call ccp4h_rule() call ccp4h_summary_beg() call ccp4h_header('Scale from Wilson Plot','wilson',3) call ccp4h_summary_end() call ccp4h_pre_beg() C C---- Calculate points for wilson plot C C Plot log(fp*fp/ff) = log (1/k) -2B*STHSQ C = log (1/k) - B*S/2 C JJ= 0 WRITE(6,'(/////,1X,a)') ' ******* Wilson Plots *******' WRITE(6,6543) 6543 FORMAT(/, +' Nref is the number of observed reflections in a',/, +' hemisphere of reciprocal space.',/, +' N_unq is an estimate of the number of possible reflections',/, +' in an assymmetric unit of reciprocal space ',/, +' ( Nref should be approximately equal to Nsymp*N_unq)',/, +' Mn(ff) is the expected value of f**2',/, +' Mn(s - resln) ', $ 'is the average value of 4(sin theta/lambda)**2 and ',/, +' the corresponding resolution limit.',/, +' Mn(fobssq) is the average value of Fobs**2') WRITE(6,6533) 6533 FORMAT(//, +' If the reflections which were not measured were all weak,', +' then Mn(Fobs**2)',/, $' is better estimated using all possible', +' reflections N_unq',/, $' (Option WILSON ALL). THIS SHOULD NOT', $' NORMALLY BE USED',/, +' ln(Mn((Fo**2)1))/Mn(ff) uses the average derived from Nref,' +,/, +' ln(Mn((Fo**2)2))/Mn(ff) uses the average derived from N_unq.') call ccp4h_graph_beg(700,300) WRITE(6,'(///,3(A,/))') $ ' $TABLE: Wilson Plot:', $ ' $GRAPHS:Wilson Plot:A:5,8,10:',' $$' WRITE(6,'(12(a,/))') 1 'Range', 1 ' N_obs(P1)', 1 ' N_possible(P1)', 1 ' Mn(ff_allatoms)', 1 ' 1/resol^2', 1 ' Resolution', 1 ' Mn(Fsquared/n_obs) ', 1 ' ln(Mn(Fsquared/n_obs)/Mn(ff) ', 1 ' Mn(Fsquared/n_poss) ', 1 ' ln(Mn(Fsquared/n_poss)/Mn(ff) $$ ', 1 ' i nref N_unq Mn(ff) 1/resol^2 Mn(resln) Mn(FF/n_obs)', 1 ' ln(Mn(FF/n_obs)/Mn(ff)) Mn(FF/n_poss)', 1 ' lnMn(FF/n_poss)/Mn(ff) ', $ ' $$' C C IOBS1 = 0 IOBS2 = NRANGE DO 301 I=1,NRANGE C DO 301 I=1,NBIN C IF(I.EQ.1)NOPT = NRFUNQ(I) C IF(I.GT.1)NOPT = NRFUNQ(I) - NRFUNQ(I-1) NOPT = NRFUNQ(I) IF(SW(I).EQ.0.0 .OR. SN(I) .EQ.0.0) GO TO 301 C JJ = JJ + 1 SW(I)= SW(I)/SN(I) SR(I) = SR(I)/SN(I) FFAV1 = SI(I)/SN(I) IF(IOBS1 .EQ. 0) IOBS1 = I IF(SR(I).LT.SMINSC) IOBS1= I + 1 IF(SR(I).LE.SMAXSC) IOBS2= I IFFAV1 = NINT(FFAV1) IF (SR(I) .GT. 0.0) RESI = 1.0/SQRT( SR(I) ) NRFF = NINT(SN(I)) IF(SW(I).GT.0.00001 .AND. FFAV1 .GT. 0.0) THEN FF1(I) = ALOG(FFAV1/SW(I)) ELSE IF(FFAV1 .LE. 0.0) THEN WRITE (6,6101) ' Negative or zero mean I, bin', I, $ ', resolution ', resi,'A', $ ' Change high resolution limit to exclude non-existent data!' 6101 FORMAT (//A,I4,A,F8.2,A,//,A//) CALL CCPERR(1, $ '*** Data beyond useful resolution limit ***') ENDIF ENDIF FFAVV1(I) = FFAV1 C C---- Now assume all missing Fobs = 0 - ie will equal C Sig(fobs**2)/ N_unq*nsymp C No need to change - they won't alter much C for each range.... C IF (NOPT .NE. 0.0) FFAV2 = SI(I) /(NOPT) IFFAV2 = NINT(FFAV2) FFAVV2(I) = FFAV2 FF2(I) = ALOG(FFAV2/SW(I)) C C---- Divide 4sinsq/lsq by 2.0 for lsq - see equation for b/scale C XYAV(1,I) = SR(I) XY(1,I) = SR(I) XLSQ(I) = XY(1,I)/2.0 SIG(I) = 1.0 C C---- First wilson plot uses ff1(...) C XYAV(2,I) = FFAVV1(I) XY(2,I) = FF1(I) YLSQ(I) = XY(2,I) WRITE(6,'(I4,2I6,F11.0,F9.4,F6.2,I10,4X,F9.5,4X,I10,4X,F9.5)') 1 I,NRFF,NOPT,SW(I),SR(I),RESI,IFFAV1,FF1(I),IFFAV2,FF2(I) 301 CONTINUE C WRITE (6,'(A)') ' $$' call ccp4h_graph_end() C WRITE (6,9876) NREFWL 9876 FORMAT (//,' A total of ',I9, + ' reflections were included in the Wilson plot',/) C TITLEX = ' 4*sinsq/lamdbasq ' TITLEY = ' mean(FPsq) ' LEGEND = ' Mean Fsq calculated using observed reflns only' IF(IPLFSQ.EQ.1)WRITE(6,'(''1'')') C C IF(IPLFSQ.EQ.1)CALL PLOTR(XYAV,NRANGE,0,'*',TITLEX,TITLEY,LEGEND) C C TITLEX = ' 4*sinsq/lamdbasq' TITLEY = ' log(FPsq/Mn(ff)) ' LEGEND = ' WILSON PLOT (observed reflections only)' WRITE(6,'(a)') '1' C C CALL PLOTR(XY,NRANGE,0,'*',TITLEX,TITLEY,LEGEND) C C---- Calculate points for wilson plot 1. C Plot log(fp*fp/ff) = log (1/k) -2B*STHSQ C = log (1/k) - B*S/2 C CALL LSQ(XLSQ,YLSQ,IOBS2,IOBS1,SIG,1,A,B,SIGA,SIGB,CHI2,Q) C C SIGRMS = -SIGB*CORRN RESI1 = 1.0/SQRT( SR(IOBS1) ) RESI2 = 1.0/SQRT( SR(IOBS2) ) call ccp4h_summary_beg() WRITE (6,'(///,a,i3,a,i3,/,a,2f8.4)') 1 ' WILSON PLOT for Ranges ',IOBS1,' - ',IOBS2, 1 ' Resolution range:',RESI1,RESI2 call ccp4h_summary_end() call ccp4h_pre_beg() WRITE (6,FMT='(A,F15.6)') ' LSQ Line Gradient = ',B WRITE (6,FMT='(A,E15.4)') 1 ' Uncertainty in Gradient = ',SIGB WRITE (6,FMT='(A,E15.4)') ' X Intercept = ',A WRITE (6,FMT='(A,E15.4)') 1 ' Uncertainty in Intercept = ',SIGA C C---- Calculate points for wilson plot C C Plot log(fp*fp/ff) = log (1/k) -2B*STHSQ C = log (1/k) - B*S/2 C call ccp4h_summary_beg() WRITE(6,'(///,a,/,a)') 1 ' For a wilson plot B = - gradient', 1 ' SCALE = exp( - intercept).' C C WILBT(1) = - B WILSC(1) = EXP(-A) WRITE (6,1120) WILBT(1),WILSC(1) 1120 FORMAT(//,' Least squares straight line gives: B =', 1 F7.3,8X,'SCALE =',F10.5, 1 /,' where F(absolute)**2 = ', 1 'SCALE*F(observed)**2*EXP(-B*2*SINTH**2/L**2) ',///) call ccp4h_summary_end() C C---- Now assume all missing Fobs = 0 - ie will equal C Sig(fobs**2)/ N_unq*nsymp C No need to change - they won't alter much C for each range.... C C---- Second wilson plot uses ff2(...) C IF (IWITCH .EQ. 2) THEN DO 401 I=1,NRANGE XYAV(2,I) = FFAVV2(I) XY(2,I) = FF2(I) YLSQ(I) = XY(2,I) 401 CONTINUE C C TITLEX = ' 4*sinsq/lamdbasq' TITLEY = ' mean(FPsq) ' LEGEND = ' Mean Fsq calculated using n= all possible reflns' IF(IPLFSQ.EQ.1)WRITE(6,'(a)') '1' IF(IPLFSQ.EQ.1) 1 CALL PLOTR(XYAV,NRANGE,0,'*',TITLEX,TITLEY,LEGEND) C C TITLEX = ' 4*sinsq/lamdbasq' TITLEY = ' log(FPsq/Mn(ff)) ' LEGEND = + ' WILSON PLOT (all possible reflections)' WRITE(6,'(a)') '1' CALL PLOTR(XY,NRANGE,0,'*',TITLEX,TITLEY,LEGEND) C C---- Calculate points for wilson plot C Plot log(fp*fp/ff) = log (1/k) -2B*STHSQ C = log (1/k) - B*S/2 C CALL LSQ(XLSQ,YLSQ,IOBS2,IOBS1,SIG,1,A,B,SIGA,SIGB,CHI2,Q) C C SIGRMS = -SIGB*CORRN WRITE (6,'(///,a)')' WILSON PLOT 2 ' WRITE (6,FMT='(A,F15.6)') ' LSQ Line Gradient = ',B WRITE (6,FMT='(A,E15.4)') 1 ' Uncertainty in Gradient = ',SIGB WRITE (6,FMT='(A,E15.4)') ' X Intercept = ',A WRITE (6,FMT='(A,E15.4)') 1 ' Uncertainty in Intercept = ',SIGA C C---- Calculate points for wilson plot C C Plot log(fp*fp/ff) = log (1/k) -2B*STHSQ C = log (1/k) - B*S/2 C WILBT(2) = - B WILSC(2) = EXP(-A) WRITE (6,1121) WILBT(2),WILSC(2) 1121 FORMAT(//,' Least squares straight line gives: B =', 1 F7.3,8X,'SCALE =',F10.5, 1 /,' Where F(Absolute)**2 = ', 1 'SCALE*F(observed)**2*EXP(-B*2*SINTH**2/L**2) ',///) C ENDIF C call ccp4h_pre_end() RETURN END C SUBROUTINE WILPRE(NSYM,CELL) C ============================ C C IMPLICIT NONE C INTEGER NSYM REAL CELL(6) C C INTEGER NBIN,NMON,IPLFSQ,NK REAL RMAX,RMIN,SMIN,SMAX,SMINSC,SMAXSC,RMINSC,RMAXSC, $ VPA COMMON /WILKEY/ NBIN,NMON,RMAX,RMIN,SMIN,SMAX,SMINSC, + SMAXSC,RMINSC,RMAXSC,VPA, + IPLFSQ,NK REAL DSS(100),SI,SN,SR,SW,XY,XYAV,FFS,PI,DELDSS,F000, + VOL,RR,CELLAS,RVOL,FRAC INTEGER NRFUNQ,NATMS,NELEC,ISMAX COMMON /PREWIL/ SI(100),SN(100),SR(100),SW(100), + NRFUNQ(100),XY(2,100),XYAV(2,100), + FFS(500),PI,DELDSS,NATMS,NELEC,F000, + ISMAX,VOL,RR(3,3,6),CELLAS(6),RVOL,FRAC C CHARACTER NW*2 CHARACTER NW*4,ATID*4 COMMON /WILCHR/ NW(20) INTEGER NC REAL AF(4),BF(4) COMMON /WILFRM/ NC(20),KatmWghts(20),FormFactAcoeff(5,20), + FormFactBcoeff(4,20) C COMMON /WILFRM/ NC(8),NO(8),AL(8),AS(8),BL(8),BS(8),CL(8) C C Locals INTEGER LUNOUT,I,J LOGICAL IFILFF REAL S0,S1,SS,FZ C LUNOUT = 6 C C C WRITE (LUNOUT,8888) 8888 FORMAT(' In Subroutine WILPRE',//) C C IF (NBIN.GT.60) THEN NBIN = 60 WRITE(LUNOUT,'(A,I5)') 1 ' Maximum number of Ranges for analysis reset to ',NBIN END IF C C CALL RBFRO1(CELL,VOL,RR) CALL RBRCEL(CELLAS,RVOL) C C---- Now check if Wilson plot information been input. If not C guess the number of residues assuming 50% solvent in the cell C and assigning a average volume of 157A**3 to each residue. C Rough but not totally unreasonable - based on 5+1.35+1.5+8 atoms per residue; C Volume per atom ~ 10A**3 C NRES = 999 IF (NK.EQ.0) THEN NRES = 0.5*VOL/(NSYM*157) C NW(1) = 'C' NC(1) = 5 * NRES NW(2) = 'N ' NC(2) = NINT(1.35 * NRES) NW(3) = 'O ' NC(3) = NINT(1.5 * NRES) NW(4) = 'H ' NC(4) = 8 * NRES NK = 4 NRES = -999 END IF C C C WRITE (LUNOUT,FMT=9010) RMIN,RMAX,SMIN,SMAX 9010 FORMAT (/, + ' Resolution limits in As = ',2F10.2,/, + ' as 4sinsq/lsq = ',2F10.5,/) WRITE (LUNOUT,FMT=9011) RMINSC,RMAXSC,SMINSC,SMAXSC 9011 FORMAT (/, + ' Resolution limits used for scaling in As = ',2F10.2,/, + ' as 4sinsq/lsq = ',2F10.5,/) call ccp4h_pre_end() C C C IFILFF = .TRUE. DO 5 IFMFC = 1,NK ATID = NW(IFMFC) C C ************************* CALL SFREAD2(ATID, + 5, + AF, + BF, + CF, + IWT, + IELEC, + CU, + MO, + Ifail) C ************************* C C ***************************************** IF (Ifail.EQ.-1) + CALL CCPERR(1,' No match for atom label(KEYIN1)') C ***************************************** C KatmWghts(IFMFC) = IWT FormFactAcoeff(1,IFMFC) = AF(1) FormFactAcoeff(2,IFMFC) = AF(2) FormFactAcoeff(3,IFMFC) = AF(3) FormFactAcoeff(4,IFMFC) = AF(4) FormFactBcoeff(1,IFMFC) = BF(1) FormFactBcoeff(2,IFMFC) = BF(2) FormFactBcoeff(3,IFMFC) = BF(3) FormFactBcoeff(4,IFMFC) = BF(4) FormFactAcoeff(5,IFMFC) = CF 5 CONTINUE C NumFormFactors = NK C C CALL ATREC(AL,AS,BL,BS,CL,NK,NC,NW) C C C IF(IPLFSQ.EQ.0)WRITE (LUNOUT,'(/,A)') ' No plot of Fsq ' C C---- Grid on 4sinsq/lamdasq C DELDSS = SMAX/NBIN C C---- How many electrons in assym unit?? C NATMS = 0 NELEC = 0 C C DO 151 I=1,NK C FFF = AL(I) + BL(I) + CL(I) C NO(I) = NINT(FFF) C NO(I) = KatmWghts(IFMFC) NELEC=NELEC+ NSYM*KatmWghts(I)*NC(I) NATMS=NATMS+NSYM*NC(I) 151 CONTINUE C C F000 = NELEC C C FRAC = REAL(NATMS)*VPA/VOL C C call ccp4h_rule() call ccp4h_header('Volume, Solvent Content etc','volume',3) call ccp4h_pre_beg() WRITE (LUNOUT,'(///,1X,A)') + ' **** Volume solvent content etc ***' IF(NRES.LT.0) WRITE (LUNOUT,'(///,1X,A)') + ' **** NOTE: Following figures ASSUME 50% solvent', + ' No Contents information given in input ****' C C + (NW(I),NC(I),NO(I),AL(I),AS(I),BL(I),BS(I),CL(I),I=1,NK) WRITE (LUNOUT,2800) + (NW(I),NC(I),KatmWghts(I),(FormFactAcoeff(J,I), + FormFactBcoeff(J,I),J=1,4),FormFactAcoeff(5,I),I=1,NK) 2800 FORMAT(/, + ' Asymmetric Unit Contents ',5X, + ' Scattering Factor Constants',/, + ' Atom number in A.U. Atomic number',5X, + '(F = AA*EXP(-A*RHO) + BB*EXP(-B*RHO) + .. + CC)',/, + (4X,A4,2I13,10X,9F9.3)) WRITE (LUNOUT,2820) + VPA,NATMS,VOL,F000,FRAC,RMIN,RMAX,DELDSS 2820 FORMAT(/, 1 ' Volume per atom = ',F13.1, + ' A**3',/, 1 ' total number of atoms in unit cell = ',I13,/, 2 ' unit cell volume = ',F13.1,/, 2 ' F000 = ',F13.1,/, 3 ' fraction of unit cell occupied by atoms = ',F13.3, $ ' =====',/, 4 ' starting resolution = ',F13.2,/, 5 ' finishing resolution = ',F13.2,/, 5 ' resolution increment for plotting = ',F13.2,//) call ccp4h_pre_end() C C---- harvest C _cell.solvent_fraction C CALL Hcell_percent_solvent(100.0 - 100.0*FRAC) C C---- harvest C C---- Number of expected reflections per shell C S0 = 0.0 S1 = DELDSS C C DO 160 I = 1,NBIN + 1 DSS(I) = (I-0.5)*DELDSS FFS(I) = 0.0 SS = DSS(I)/4.0 C C DO 156 J=1,NK FZ = FormFactAcoeff(5,J) DO 157 K=1,4 FZ = FZ + FormFactAcoeff(K,J)*EXP(-FormFactBcoeff(K,J)*SS) C FZ = (AL(J)*EXP(-AS(J)*SS) + BL(J)*EXP(-BS(J)*SS) + CL(J))**2 157 CONTINUE FFS(I) = FFS(I) + FZ*FZ*REAL(NSYM*NC(J)) 156 CONTINUE C NRFUNQ(I) = NINT(0.6667*PI*DSS(I)**1.5 /( RVOL*NSYM)) C WRITE(LUNOUT,1667) C + S0,S1,DSS(I),FFS(I),NRFUNQ(I) C 1667 FORMAT( F9.4,'-',F6.4, F8.4, F16.2,I10) C S0 = S1 C S1 = S0 + DELDSS 160 CONTINUE RETURN END C C-------------------------------- SUBROUTINE SCNDER(X,Y,N,Y2) C-------------------------------- C Returns second derivatices (Y2) of the interpolating function for y(x) C whose N values Y are tabulated for points X. C INTEGER NMAX,N,K,I PARAMETER (NMAX=100) REAL X(N),Y(N),Y2(N),U(NMAX),UN,QN,P,SIG C Y2(1)=0. U(1)=0. C DO 11 I=2,N-1 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2. Y2(I)=(SIG-1.)/P U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 11 CONTINUE C QN=0. UN=0. Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.) C DO 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE C RETURN END C--------------------------------------- SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y) C--------------------------------------- C Returns a cubic-spline interpolated value of function y(x) C whose N values YA are tabulated for points XA. C The second derivatives of the interpolating function are input in Y2A. C INTEGER N,K,KHI,KLO REAL XA(N),YA(N),Y2A(N),H,A,B,Y,X C KLO=1 KHI=N C 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF(XA(K).GT.X)THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF C H=XA(KHI)-XA(KLO) IF (H.EQ.0.0) CALL CCPERR(1,'SPLINT: BAD XA INPUT.') A=(XA(KHI)-X)/H B=(X-XA(KLO))/H Y=A*YA(KLO)+B*YA(KHI)+ * ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6. RETURN END C The following routines (FALLOFF,PLOTB,PTHICK,TRICART) kindly C donated by Yorgo Modis, from his program FALLOFF. Incorporated C into TRUNCATE by Martyn Winn. SUBROUTINE FALLOFF c----This program calculates the fall-off of rms F-values as a function of c----(sinth/lab)**2. This is done in 3 orthogonal directions. An overall c----fall-off of all reflections is also calculated. c---- C MDW: calculation moved to main section. This writes out tables, C and calls PLOTB to make plot. IMPLICIT NONE INTEGER KNBINS PARAMETER (KNBINS=60) INTEGER ncols PARAMETER (ncols=19) INTEGER MAXHIS PARAMETER (MAXHIS=30) integer i,j integer ihsum,ksum,lsum,osum,icheck real xx logical KYFALL,plotx,plotfall CHARACTER CTPRGO*1,CTPRGI*1,LSPRGO*30,LSPRGI*30,HISOUT*80 CHARACTER DATEMT*8,VERSNX*10,TITIN*70,TITOUT*70,LBOLIN*400 real somdir(3,KNBINS),somov(KNBINS),somsddir(3,KNBINS), + somsdov(KNBINS),sthstap,cone,somtot(4) integer numdir(3,KNBINS),numov(KNBINS),nmax,nmin,numtot(4) COMMON /FALL/KYFALL,PLOTX,PLOTFALL,cone COMMON /COORD/nmin,nmax,sthstap,somdir,somov,somsddir,somsdov, + numdir,numov COMMON /CHRMTZ/LSPRGI(16),CTPRGI(16), + LSPRGO(ncols),CTPRGO(ncols),HISOUT(MAXHIS), + VERSNX,TITIN,TITOUT,DATEMT,LBOLIN write(6,'(/,A,A,/,A,A)') + ' ANALYSIS OF THE ANISOTROPY OF THE DATA ', + 'ACCORDING TO THE FALLOFF PROCEDURE.', + ' ---------------------------------------', + '-----------------------------------' write(6,'(/)') write(6,*)'Direction 1 is perpendicular to b* and Direction 3' write(6,*)'Direction 2 is along b*' write(6,*)'Direction 3 is perpendicular to a* and b*' write(6,'(/)') C C Maire altered format statement so < 132 characters C C--- xloggraph call ccp4h_graph_beg(0,0) WRITE(6,6400) TITOUT(1:40) 6400 FORMAT(/1X, $' $TABLE: Anisotropy analysis (FALLOFF). ',A,':',/, $' $GRAPHS:Mn(F) v. resolution:A:1,2,3,4,5:', $': Mn(F/sd) v. resolution:A:1,6,7,8,9:', $': No. reflections v. resolution:A:1,10,11,12,13:',/,' $$',/, $' 1/resol^2 Mn(F(d1)) Mn(F(d2))', $' Mn(F(d3)) Mn(F(ov)) Mn(F/sd(d1)) Mn(F/sd(d2))', $' Mn(F/sd(d3)) Mn(F/sd(ov)) N(d1) N(d2) ', $' N(d3) N(ov)',/,' $$',/,' $$') somtot(1)=0 somtot(2)=0 somtot(3)=0 somtot(4)=0 numtot(1)=0 numtot(2)=0 numtot(3)=0 numtot(4)=0 do i=nmin,nmax do j=1,3 if(numdir(j,i).eq.0) then somdir(j,i)=0 somsddir(j,i)=0 else somtot(j)= somtot(j) + somdir(j,i) numtot(j)= numtot(j) + numdir(j,i) somdir(j,i)=somdir(j,i)/numdir(j,i) somsddir(j,i)=somsddir(j,i)/numdir(j,i) endif enddo if(numov(i).eq.0) then somov(i)=0 somsdov(i)=0 else somtot(4)= somtot(4) + somov(i) numtot(4)= numtot(4) + numov(i) somov(i)=somov(i)/numov(i) somsdov(i)=somsdov(i)/numov(i) endif write(6,15) i*(4*sthstap),somdir(1,i),somdir(2,i), + somdir(3,i),somov(i),somsddir(1,i),somsddir(2,i), + somsddir(3,i),somsdov(i),numdir(1,i),numdir(2,i), + numdir(3,i),numov(i) 199 continue enddo 15 format(1X,f7.5,8(1x,f8.2),4(1x,i6,2X)) WRITE(6,'(A)') ' $$' call ccp4h_graph_end() icheck = 0 somtot(4) = somtot(4)/numtot(4) do j=1,3 if(numtot(j) .gt. 0) then somtot(j) = somtot(j)/numtot(j) if(abs(somtot(j)/somtot(4) - 1.0) .gt. 0.2) icheck = 1 endif end do write(6,'(/,A,3F9.2,4x,F9.2)') + ' Average F (d1 d2 d3) + overall average:',somtot c if(icheck .eq. 1) call ccperr(2, +' ***Beware-serious ANISOTROPY; data analyses may be invalid ***') ihsum=0 ksum=0 lsum=0 osum=0 do 1075 i=nmin,nmax ihsum=ihsum+numdir(1,i) ksum=ksum+numdir(2,i) lsum=lsum+numdir(3,i) osum=osum+numov(i) 1075 continue write(6,6010) 1 cone,ihsum,cone,ksum,cone,lsum,osum 6010 format (//' number A-AX reflections less than ',f4.1, 1' degrees from dir1',i10, 1 /' number B-AX reflections less than ',f4.1, 1' degrees from dir2',i10, 1 /' number C-AX reflections less than ',f4.1, 1' degrees from dir3',i10, 1 /' number overall reflections ',i10) IF (PLOTFALL) THEN write(6,*)'Direction 1 is plotted as a thick line' write(6,*)'Direction 2 is plotted as a hollow line with boxes at +regular intervals of resolution' write(6,*)'Direction 3 is plotted as a thin line' write(6,'(/)') C ***** CALL PLOTB C ***** C XX = 10.0 C C *************************************** CALL GSMVTO(XX,0.0) CALL GSENDP CALL GSSTOP C *************************************** ENDIF write(6,'(/)') return end C C ================ SUBROUTINE PLOTB C ================ C INTEGER KNBINS PARAMETER (KNBINS=60) INTEGER ncols PARAMETER (ncols=19) INTEGER MAXHIS PARAMETER (MAXHIS=30) C .. C .. Scalars in Common .. integer nmax,nmin,NRANGE logical KYFALL,plotx,PLOTFALL real x_scale,y_scale,RANGE,cone C C .. C .. Local Scalars .. REAL CSIZE,DBMULT,SDD,DOTSMM,HCHARB,HSYMB, + PAPLEN,PAPWID,PI,PI2,RDELTAY,TSIZE,TSTEP,X0,XT,XX,XYZ_TIC, + Y0,YT,YY INTEGER I,JDO,KBSTEP,KCOUNT,KRSTEP,MIXCLR,NAMRB,NFONT, + NLINES,NTHICK CHARACTER CBWORK*5,reslab*6 C .. C .. Local Arrays .. REAL BOX_X(5),BOX_Y(5) C .. C .. Arrays in Common .. CHARACTER CTPRGO*1,CTPRGI*1,LSPRGO*30,LSPRGI*30,HISOUT*80 CHARACTER DATEMT*8,VERSNX*10,TITIN*70,TITOUT*70,LBOLIN*400 real somdir(3,KNBINS),somov(KNBINS),somsddir(3,KNBINS), + somsdov(KNBINS) integer numdir(3,KNBINS),numov(KNBINS) C .. C .. External Subroutines .. EXTERNAL CCPERR,GSANCU,GSBSIZ,GSCENC,GSCETX,GSCOLR,GSCSPU,GSDTRN, + GSDWTO,GSFONT,GSINIT,GSLNWT,GSMIXC,GSMVTO,GSORGD,GSPICT, + GSSCLC,GSUROT,PTHICK C .. C .. Common blocks .. COMMON /KEYRSC/RANGE,NRANGE COMMON /FALL/KYFALL,PLOTX,PLOTFALL,cone COMMON /COORD/nmin,nmax,sthstap,somdir,somov,somsddir,somsdov, + numdir,numov COMMON /CHRMTZ/LSPRGI(16),CTPRGI(16), + LSPRGO(ncols),CTPRGO(ncols),HISOUT(MAXHIS), + VERSNX,TITIN,TITOUT,DATEMT,LBOLIN c C .. Data statements .. C DATA BOX_X/20.0,270.0,270.0,20.0,20.0/ DATA BOX_Y/10.0,10.0,140.0,140.0,10.0/ C .. C C X_SCALE = 0.5 Y_SCALE = 0.5 FMAX = 0.0 C MIXCLR = 1 PI = ATAN2(1.0,1.0)*4.0 PI2 = PI/2.0 PAPWID = 328.0*X_SCALE PAPLEN = 508.0*Y_SCALE NFONT = 1 DOTSMM = 10.0 NTHICK = 1 C C---- Centered symbol height ('mms) B-Plot C HSYMB = 1.5 C C---- Character height (mms) B-Plot C HCHARB = 3.5 CALL GSINIT('PLOT') CALL GSMIXC(MIXCLR) CALL GSBSIZ((PAPWID+1.0), (PAPLEN+1.0)) CALL GSDTRN(DOTSMM,DOTSMM) CALL GSPICT CALL GSFONT(NFONT) C C---- Centred characters with uniform spacing C CALL GSCENC(1) CALL GSCSPU(1) C C---- Plot X C IF (PLOTX) THEN CALL GSORGD(0.0,0.0) CALL GSUROT(0.0,PI2) ELSE C C---- Plot Y C CALL GSORGD(PAPWID,0.0) CALL GSUROT(PI2,PI) END IF C C CSIZE = HCHARB*X_SCALE CALL GSSCLC(CSIZE,CSIZE) CALL GSLNWT(NTHICK) CALL GSCOLR(1) C C---- Draw Box C CALL GSMVTO(BOX_X(1)*X_SCALE,BOX_Y(1)*Y_SCALE) C C DO 10 JDO = 2,5 CALL GSDWTO(BOX_X(JDO)*X_SCALE,BOX_Y(JDO)*Y_SCALE) 10 CONTINUE C C---- Title C CALL GSLNWT(2) CALL GSCOLR(5) TSIZE = HCHARB*1.60*X_SCALE CALL GSSCLC(TSIZE,TSIZE) CALL GSMVTO(143.0*X_SCALE,145.0*Y_SCALE) CALL GSANCU(143.0*X_SCALE,145.0*Y_SCALE) CALL GSCETX(TITOUT,2) C C---- Find max C DO 20 I = nmin,nmax do j=1,3 FMAX = MAX(FMAX,somdir(j,I)) enddo 20 CONTINUE C C NFONT = 1 NTHICK = 1 CALL GSLNWT(NTHICK) CALL GSCOLR(1) C C C CSIZE = 2.75*X_SCALE CALL GSMVTO(16.0*X_SCALE,15.0*Y_SCALE) CALL GSDWTO(20.0*X_SCALE,15.0*Y_SCALE) C C C---- first C k=NRANGE-1 RDELTAY = 240.0/k C C---- loop res labels C KRSTEP = (nmax-nmin)/NRANGE KCOUNT = 1 C C---- zero C CALL GSMVTO(15.0*X_SCALE,15.0*Y_SCALE) CALL GSANCU(15.0*X_SCALE,15.0*Y_SCALE) CSIZE = HCHARB*X_SCALE CALL GSSCLC(CSIZE,CSIZE) CALL GSCETX(' 0 ',3) C C---- C FT = MOD(FMAX,10.0) FMAXT = FMAX + FT + 0.51 NAMRB = NINT(FMAXT/10.0)*10 DBMULT = 120.0/NAMRB KBSTEP = FMAX/10 c c DO 500 JDO = KBSTEP,NAMRB,KBSTEP CBWORK = ' ' IF (JDO.LE.9) THEN WRITE (CBWORK,FMT=6000) JDO 6000 FORMAT (' ',I1,' ') GO TO 40 ELSE IF (JDO.LE.99) THEN WRITE (CBWORK,FMT=6002) JDO 6002 FORMAT (' ',I2,' ') GO TO 40 else if (jdo.le.999) then write (cbwork,fmt=6003) jdo 6003 format (' ',i3,' ') go to 40 else if (jdo.le.9999) then write (cbwork,fmt=6004) jdo 6004 format (i4,' ') go to 40 else if (jdo.le.99999) then write (cbwork,fmt=6005) jdo 6005 format (i5) csize=csize/1.3 CALL GSSCLC(CSIZE,CSIZE) go to 40 endif 40 CONTINUE C C XYZ_TIC = (DBMULT*JDO) + 15.0 CALL GSMVTO(15.0*X_SCALE,XYZ_TIC*Y_SCALE) CALL GSANCU(15.0*X_SCALE,XYZ_TIC*Y_SCALE) CALL GSCETX(CBWORK,3) CALL GSMVTO(16.0*X_SCALE,XYZ_TIC*Y_SCALE) CALL GSDWTO(20.0*X_SCALE,XYZ_TIC*Y_SCALE) csize= HCHARB*X_SCALE CALL GSSCLC(CSIZE,CSIZE) 500 CONTINUE call gsmvto(18.*x_scale,(dbmult*(namrb+kbstep)+15.)*y_scale) call gsancu(18.*x_scale,(dbmult*(namrb+kbstep)+15.)*y_scale) call gscetx('Mn(|F|)',3) DO 501 sDO = nmin*(4*sthstap),nmax*(4*sthstap), +16*sthstap sdd=sqrt(1/sdo) write (reslab,fmt=6010) sdd 6010 format(f6.2) dsmult=231.4/(nmax*(4*sthstap)) XYZ_TIC =(DsMULT*sDO) + 35.0 CALL GSMVTO(XYZ_TIC*X_SCALE,2.0*Y_SCALE) CALL GSANCU(XYZ_TIC*X_SCALE,2.0*Y_SCALE) CALL GSCETX(reslab,2) CALL GSMVTO(XYZ_TIC*X_SCALE,7.0*Y_SCALE) CALL GSDWTO(XYZ_TIC*X_SCALE,10.0*Y_SCALE) c sdp=(nmax+6)*(4*sthstap) c call gsmvto((1.*dsmult*sdp+0.)*x_scale,3.0*y_scale) c call gsancu((1.*dsmult*sdp+0.)*x_scale,3.0*y_scale) c call gscetx('d (Angstrom)',3) 501 continue call gsmvto(xyz_tic*x_scale+5.,3.*y_scale) call gsancu(xyz_tic*x_scale+5.,3.*y_scale) call gscetx('d (Angstrom)',1) C c-- fs C CALL GSLNWT(3) CALL GSCOLR(4) C do j=1,3 NLINES = 8 TSTEP = 0.075 Y0 = ((somdir(1,1)*DBMULT+15.0)*Y_SCALE) X0 = (25.0*X_SCALE) c c DO 90 i = nmin,nmax-1 IF (somdir(j,i).EQ.-999.0) GO TO 90 if (somdir(j,i).eq.0.0) go to 90 XX = (((i-1)*RDELTAY+25.0)*X_SCALE) if (i.eq.nmin) then Y0 = ((somdir(j,1)*DBMULT+15.0)*Y_SCALE) x0=(25.0*X_SCALE) endif YY = ((somdir(j,i)*DBMULT+15.0)*Y_SCALE) XT = XX YT = YY c c-make line for dir2&3 different thickness c if (j.eq.2) then call pthick(x0,y0,xx,yy,0.75,2) call pthick(x0-0.53,y0-0.5,x0+0.53,y0-0.5,tstep,4) call pthick(x0-0.5,y0-0.5,x0-0.5,y0+0.5,tstep,4) call pthick(x0-0.53,y0+0.5,x0+0.53,y0+0.5,tstep,4) call pthick(x0+0.5,y0-0.5,x0+0.5,y0+0.5,tstep,4) else if (j.eq.3) then call pthick(x0,y0,xx,yy,tstep,2) else CALL PTHICK(X0,Y0,XX,YY,TSTEP,NLINES) endif X0 = XT Y0 = YT 90 CONTINUE enddo C C RETURN END C C =========================================== SUBROUTINE PTHICK(X0,Y0,X1,Y1,TSTEP,NLINES) C =========================================== C C Draw line from (X0,Y0) to (X1,Y1) C C Parameters to control line thickness as function of intensity C TSTEP intensity step for each increase in thickness C NLINES maximum number of lines C TSTEP offset step in plotter units C C C C C---- always draw central line C C .. Scalar Arguments .. REAL TSTEP,X0,X1,Y0,Y1 INTEGER NLINES C .. C .. Local Scalars .. REAL STEP,X,X2,XM,XT,Y,Y2,YM,YT INTEGER JDO,NOFF LOGICAL HORIZ C .. C .. External Subroutines .. EXTERNAL GSDWTO,GSMVTO C .. C .. Intrinsic Functions .. INTRINSIC ABS,MOD C .. CALL GSMVTO(X0,Y0) CALL GSDWTO(X1,Y1) C C---- set up current point C XM = X0 YM = Y0 X = X1 Y = Y1 C C---- draw thickening lines C If line nearer horizontal, offsets will be on y and vv C HORIZ = .FALSE. C C IF (ABS(X1-X0).GT.ABS(Y1-Y0)) HORIZ = .TRUE. C C DO 10 JDO = 2,NLINES C C---- Offset alternates between + & -, increasing in size C NOFF = JDO/2 IF (MOD(JDO,2).NE.0) NOFF = -NOFF C C---- Apply offsets C STEP = NOFF*TSTEP C C IF (HORIZ) THEN C C---- . . to y C XT = XM YT = YM + STEP X2 = X Y2 = Y + STEP ELSE C C---- . . to x C XT = X0 + STEP YT = YM X2 = X + STEP Y2 = Y END IF C C---- Plot alternate lines backwards C IF (NOFF.LT.0) THEN CALL GSMVTO(XT,YT) CALL GSDWTO(X2,Y2) ELSE CALL GSMVTO(X2,Y2) CALL GSDWTO(XT,YT) END IF C C 10 CONTINUE C C---- Reset current point C XM = X YM = Y C RETURN END subroutine tricart(a1,a2,a3,alpha,beta,gamma,transf) c Calculates the matrix that transforms coordinates relative to the c triclinic axes a1, a2 and a3 to a Cartesian set of axes. a2(cart) c runs along a2, a1(cart) lies in the plane of a1 and a2 and a3(cart) c runs along a1 x a2. c -- -- C C I.e. X // b* x c, Y // b*, Z // a* x b* // c C This does not agree with any of the standard orthogonalisation C codes, and so we cannot use library functions to get transf. c c Alpha, beta, gamma must be given in radians. c Lit.: M.G.Rossman & D.M.Blow, Acta Cryst.(1962) Vol.15,24 c formula (9). c dimension transf(3,3) c1=cos(alpha) c2=cos(beta) c3=cos(gamma) s1=sin(alpha) s3=sin(gamma) cw=(c2-c1*c3)/(s1*s3) sw=sin(acos(cw)) transf(1,1)=a1*s3*sw transf(1,2)=0.0 transf(1,3)=0.0 transf(2,1)=a1*c3 transf(2,2)=a2 transf(2,3)=a3*c1 transf(3,1)=a1*s3*cw transf(3,2)=0.0 transf(3,3)=a3*s1 return end