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 Linked with new libraries and checked 26th Sept 1990 EJD. c Alterations to PARSER calls needed parameter(nsize=60000) character types (7) integer lengs (7) integer ccpe2i external search, ccpalc, ccpe2i data types /'i', 6*'r'/ CCC dimension igrid(nsize),sigfcc(nsize,5),sigfc(nsize,5), CCC 1 siGFco(nsize,5),delf(nsize,5),delft(nsize),den(nsize) isize = ccpe2i ('RSEARCH_NSIZE', nsize) C for IGRID lengs (1) = isize C SIGFCC, SIGFC, SIGFCO, DELF do 10 i=2,5 lengs (i) = 5*isize 10 continue C DELFT, DEN lengs (6) = isize lengs (7) = isize call ccpalc (search, 7, types, lengs) call ccperr (0, 'Normal termination') end C-- SEARCH PCZ.SEARCH.FORT 26/04/85 PAM C C PROGRAM 'SEARCH' C ================ C C R-FACTOR SEARCH PROGRAM WRITTEN BY ELEANOR DODSON C C FORTRAN 77 VERSION J. CAMPBELL APRIL 1985 C C SUBROUTINE SEARCH(NSIZE, IGRID, ifcc, SIGFCC, ifc, SIGFC, ifco, + SIGFCO, idelf, DELF, idelft, DELFT, iden, DEN) C all the array dimensions set above by ccpalc are dummy except for C NSIZE DIMENSION DELTSN(5),DEN(NSIZE) DIMENSION IGRID(NSIZE),NFC(5),SIGFO(5),sigfoo(5), 1 sigfc(nsize,5), 1 sigfcc(NSIZE,5),sigfco(NSIZE,5),DELF(NSIZE,5),DELFT(NSIZE) DIMENSION AH(12),AK(12),AL(12),CELLMTZ(6), 1 IN(3),ACPART(12),BCPART(12), 1 JKP(500),FKP(500) DIMENSION IXFPMN(30),IXFPMX(30),IXFCMN(30) DIMENSION SINA(10241),COSA(8193) c Arrays for reflection selection c ie rzone -h +k +l = 3n needs irzone(j,1:5) = -1 1 1 3 0 c ie rzone l = 2n +1 needs irzone(j,1:5) = 0 0 1 2 1 dimension irzone(10,5) COMMON /ATSYM/rsym1(4,4,192),Rsym2(4,4,192),NSYM,PERM(4,4), 1 MATOMS,LATOMS c put this in to transfer nsymch no to print sr c COMMON/io/ijnk1,ijnk2,ijnk3,nspg0,nsymch,ijnk4,nsecax COMMON /IOFFT/ injnk0,ijnk1,ijnk2,NSYMCH,ijnk3,NSECAX COMMON /mphdr/ NXYZ(3),IXMIN,IXMAX,IYMIN,IYMAX,IZMIN,IZMAX, . NUSED6,NUSED7,NUSED9,NUSED10,NUSED11,IDUMMY(3),NUSED12 . ,CELLm(6),lspGRP,RHMIN,RHMAX,Rhmean,rhrms COMMON /FFTSPG/NSPGRP CHARACTER TITLE*80,namspg*10,pgnam1*10,PGNAM2*10 COMMON/mphdrr/TITLE COMMON /MAPFMT/ ITYPe, NLINES,NCOL,Vol, F000, RHOMAX,rhomin COMMON /FFTPRM/ KP1,KP2,KP3 c COMMON/sample/NX,NY,NZ COMMON/FOUT/jnk1,jnk2,ISP,IDIG,ILIN CHARACTER*1 LX,LY,LN C lcf stuff PARAMETER(MCOLS=500) CHARACTER CTPRGI(MCOLS)*1,LSPRGI(MCOLS)*30 C .. INTEGER LOOKUP(MCOLS),IPRINT REAL ADATA(MCOLS) LOGICAL LOGMSS(MCOLS) C .. Parameters .. INTEGER NKEYS,NPAR PARAMETER (NKEYS=20,NPAR=200) C .. CHARACTER LINE*600 C .. CHARACTER KEY*4,CVALUE(NPAR)*4 C .. LOGICAL LEND,EOF,LPRINT REAL FVALUE(NPAR) INTEGER IBEG(NPAR),IDEC(NPAR),IEND(NPAR),ITYP(NPAR) CHARACTER KEYWRD(NKEYS)*4 C .. C C CHARACTER*1 key1 dimension mp(3) CHARACTER*1 kxyz(3),LTYPE EQUIVALENCE (SINA(2049),COSA(1)) C EQUIVALENCE (DEN(1),sigfcc(1,1)) DATA KEYwrd/'TITL','FOLI','FCLI','SCAL','RESO','BINM','LPMA', . 'AXIS','SYMM','BGRV','CENT','REST','XGRI','YGRI', . 'ZGRI','CORR','END ','LABI','RZON','NRFA'/ DATA KXYZ/'X','Y','Z'/ DATA LX,LY,LN/'X','Y','N'/ C---- NLPRGI = number of input labels C DATA NLPRGI,LSPRGI /31,'H','K','L','FP','SIGFP', + 'FC1','PC1','FC2' ,'PC2' ,'FC3' ,'PC3' ,'FC4' ,'PC4', + 'FC5','PC5','FC6' ,'PC6' ,'FC7' ,'PC7' ,'FC8' ,'PC8', + 'FC9','PC9','FC10','PC10','FC11','PC11','FC12','PC12', + 'FC', 'PHIC',469*' '/ DATA CTPRGI /'H','H','H','F','Q', + 'F','P','F' ,'P' ,'F' ,'P' ,'F' ,'P', + 'F','P','F' ,'P' ,'F' ,'P' ,'F' ,'P', + 'F','P','F','P','F','P','F','P', + 'F', 'P',469*' '/ C C C---- This code signs which input columns are essential (LOOKUP) C DATA LOOKUP/-1,-1,-1, -1,-1, -1,-1,-1,-1, 491*0/ DATA BGRVAL/0.0/, SIGFOO/5*0.0/, IBCD /0/ CALL CCPFYP CALL MTZINI nspg0=1 IFAIL=0 CALL CCPOPN(5,'DATA',5,1,LDUM,IFAIL) CALL CCPOPN(6,'PRINTER',6,1,LDUM,IFAIL) CALL CCPRCS(6,'RSEARCH','$Date: 2002/08/06 14:05:52 $') DO 995 I=4,6 CELLM(I)=90 995 CONTINUE lspgrp=1 nspgrp=1 JDEN=NSIZE IGRDLM=NSIZE ITWO10=1024 ITWO20=1024*1024 C C Set defaults TITLE=' Rsearch run ' RESMIN=10000. RESMAX=2. map=0 nrfact=15 scale=1 bfact=0 icen=0 icorr=0 IRESTT=0 NREF=0 irz=0 C NERR=0 IXTL=0 MTZIN = 1 C isysw=6 NSYM = 0 211 CONTINUE LINE = ' ' KEY = ' ' NTOK = NPAR C 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 2100 C C---- Loop over possible key-words C DO 2000 I = 1,NKEYS IF (KEY.EQ.KEYWRD(I)) GO TO 2001 2000 CONTINUE CALL CCPERR(1, ' *** UNRECOGNISED KEY WORD ***') C C 2001 continue GO TO . (2010,2020,2030,2040,2050,2060,2070, . 2080,2090,2100,2110,2120,2130,2140, . 2150,2160,2170,2175,2190,2200),i C C TITLe read title C 2010 CONTINUE IF(NTOK.GT.1)TITLE=LINE(IBEG(2):IEND(NTOK)) WRITE(6,1001)TITLE go to 211 C C C FOLIMITS , set limits on FOBS for inclusion in R factor. C 2020 continue c 2021 CONTINUE CALL GTPREA(2,FMIN,NTOK,ITYP,FVALUE) CALL GTPREA(3,FMAX,NTOK,ITYP,FVALUE) FMIN=FMIN FMAX=FMAX if(fmin.gt.fmax)then atmp=fmin fmin=fmax fmax=atmp endif go to 211 c C FCLIMITS , set limits on FCALC for inclusion in R factor. c In low symmetry space grooups (ie P21) you should exclude all c terms where one Fcpart is very small. R factor will be independent of c translation. C 2030 continue c CALL GTPREA(2,FMAXC,NTOK,ITYP,FVALUE) FMAXC=FMAXC go to 211 c C SCAL FP , Scale and B value for FP C 2040 continue IF(LINE(IBEG(2):IEND(2)) .EQ. 'FP')THEN c CALL GTPREA(3,scale,NTOK,ITYP,FVALUE) CALL GTPREA(4,BFACT,NTOK,ITYP,FVALUE) go to 211 END IF C RESOlution read resolution limits in A, if only one treat as high C resolution limit C C RESOLUTION C 2050 CONTINUE C read resolution limits in A, if only one treat as high C resolution limit ITOK = 2 C ******************************************************** CALL RDRESO(ITOK,ITYP,FVALUE,NTOK,RESMIN,RESMAX,SMIN,SMAX) C ******************************************************** GOTO 211 C C BINMAP , set Binary map code C 2060 continue c ibin=2 c nsecax controls PRNTYZ output - always use 2 nsecax=2 IF(IBIN.GT.0)IBIN=2 go to 211 C C LPMAP , set Binary map code C 2070 continue c ibcd=6 c nsecax controls map output - default is 2 nsecax=2 go to 211 C C AXIS , set sectioning axes C 2080 continue c C AXIS ORIENTATIONS (Z,X,Y) C 2081 DO 2082 I=2,4 DO 2083 J=1,3 IF(LINE(IBEG(I):IEND(I)) .EQ. KXYZ(J))THEN MP(I-1)=J go to 2085 ENDIF 2083 CONTINUE WRITE(6,2084)LINE(IBEG(2):IEND(NTOK)) 2084 FORMAT(' *** THERE IS AN ERROR IN "AXIS"',A10,/) CALL CCPERR(1,' Stopping') 2085 CONTINUE 2082 CONTINUE c map = mp(3) go to 211 C C READ SYMMETRY CARDS C SYMMETRY C 2090 CONTINUE J = 2 IF (J .GT. NTOK) THEN WRITE (6, *) ' No symmetry data !!!' GO TO 211 END IF c C SYMMetry read crystal number or space group C C ************************************************************ CALL RDSYMM(2,LINE,IBEG,IEND,ITYP,FVALUE,NTOK,NAMSPG, + NSPGRP,PGNAM2,NSYM,NSYMP,RSYM2) C ************************************************************ write(6,*)nspgrp,nsym,nsymp write(6,*)namsPG,PGNAM2 C C RSEARCH only needs primitive symmetry operators c nsym = nsymp lspgrp = nspgrp c C SET LOOKUP FLAGS FOR HKLIN WANT FC1 AC1 ... FC(NSYMP) AC(NSYMP) MM=NSYM DO 2094 I=1,5+2*nsymp LOOKUP(i)=-1 2094 CONTINUE go to 211 C BGRV, Read input background R value. c I can't remember why this might be useful Ejd C 2100 continue c 2101 CONTINUE CALL GTPREA(2,bgrval,NTOK,ITYP,FVALUE) go to 211 C Centric refln defn. c If this is set ONLY centric reflns are included in the R value calcn. C 2110 continue c ICEN=1 go to 211 c C C RESTART H K L FOR FIRST REFLECTION TO INCLUDE C 2120 CONTINUE CALL GTPint(2,nref,NTOK,ITYP,FVALUE) CALL GTPint(3,ihh,NTOK,ITYP,FVALUE) CALL GTPint(4,ikk,NTOK,ITYP,FVALUE) CALL GTPint(5,ill,NTOK,ITYP,FVALUE) IRESTT=1 go to 211 C c Grid specifications C C EXTREMELY COMPLICATED TO ALLOW GREAT FLEXIBILITY IN DESCRIBING SEARCH VOLUME. C CODES ALLOW FOR BOXES, Prisms OR PLanes C Don't type the " please!! "x" means type the CHARACTER x c C 1) Simple rectangular box c ixmin ixmax nx c iymin iymax ny c izmin izmax nz c c C 2) Prism c ixmin ixmax nx c "xmin" iymax ny OR iymin "xmax" ny OR "xmin" "xmax" ny c izmin izmax nz c c C 3) Another sort of prism c ixmin ixmax nx c "xmin" iymax ny OR iymin "xmax" ny OR "x" "x" ny c c izmin izmax nz c OR "xmin" izmax nz OR izmin "xmax" nz OR "x" "x" nz c OR "ymin" izmax nz OR izmin "ymax" nz OR "y" "y" nz c c c c Xgrid ixmin ixmax nx c 2130 continue C CALL GTPint(2,ixmin,NTOK,ITYP,FVALUE) CALL GTPint(3,ixmax,NTOK,ITYP,FVALUE) CALL GTPint(4,nx,NTOK,ITYP,FVALUE) nxyz(1)=nx go to 211 c c Ygrid may be iymin iymax ny c or "xmin" iymax ny OR iymin "xmax" ny OR "x" "x" ny c c c Icodey 0 c OR 1 OR 2 Or 3 c 2140 continue C C ITYP(I) =0 null field C =1 character string C =2 number if(ityp(4).ne.2)go to 99 CALL GTPint(4,ny,NTOK,ITYP,FVALUE) nxyz(2)=ny c if(ityp(2).eq.2 .and.ityp(3).eq.2) then CALL GTPint(2,iymin,NTOK,ITYP,FVALUE) CALL GTPint(3,iymax,NTOK,ITYP,FVALUE) ICODEY=0 endif c C if(ityp(2).eq.1 .and.ityp(3).eq.2) then key=LINE(IBEG(2):iend(2)) CALL CCPUPC(KEY) if(key.ne.'XMIN') go to 99 CALL GTPint(3,iymax,NTOK,ITYP,FVALUE) ICODEY=1 endif c C if(ityp(2).eq.2 .and.ityp(3).eq.1) then key=LINE(IBEG(3):iend(3)) CALL CCPUPC(KEY) if(key.ne.'XMAX') go to 99 CALL GTPint(2,iymin,NTOK,ITYP,FVALUE) ICODEY=2 endif C if(ityp(2).eq.1 .and.ityp(3).eq.1) then key1=LINE(IBEG(2):(ibeg(2)+1)) CALL CCPUPC(KEY1) if(key1.ne.'X') go to 99 key=LINE(IBEG(3):(ibeg(3)+1)) CALL CCPUPC(KEY) if(key1.ne.'X') go to 99 ICODEY=3 endif go to 211 c c c Zgrid may be izmin izmax nz c OR "xmin" izmax nz OR izmin "xmax" nz OR "x" "x" nz c OR "ymin" izmax nz OR izmin "ymax" nz OR "y" "y" nz c c Icodez 0 c OR 4 OR 6 Or 8 c OR 5 OR 7 Or 9 c c 2150 continue C if(ityp(4).ne.2)go to 99 CALL GTPint(4,nz,NTOK,ITYP,FVALUE) nxyz(3)=nz c if(ityp(2).eq.2 .and.ityp(3).eq.2) then CALL GTPint(2,izmin,NTOK,ITYP,FVALUE) CALL GTPint(3,izmax,NTOK,ITYP,FVALUE) ICODEZ=0 endif c C if(ityp(2).eq.1 .and.ityp(3).eq.2) then CALL GTPint(3,izmax,NTOK,ITYP,FVALUE) key=LINE(IBEG(2):iend(2)) CALL CCPUPC(KEY) icodez=100 if(key.eq.'XMIN')ICODEZ=4 if(key.eq.'YMIN')ICODEZ=5 if(icodez.eq.100)go to 99 endif c C if(ityp(2).eq.2 .and.ityp(3).eq.1) then CALL GTPint(2,izmin,NTOK,ITYP,FVALUE) key=LINE(IBEG(3):iend(3)) CALL CCPUPC(KEY) icodez=100 if(key.eq.'XMAX')ICODEZ=6 if(key.eq.'YMAX')ICODEZ=7 if(icodez.eq.100)go to 99 endif C if(ityp(2).eq.1 .and.ityp(3).eq.1) then key1=LINE(IBEG(2):(ibeg(2)+1)) CALL CCPUPC(KEY1) icodez=100 if(key1.eq.'X')icodez=8 if(key1.eq.'Y')icodez=9 if(icodez.eq.100)go to 99 endif go to 211 c C CORRelation map plotted C 2160 continue Icorr=1 go to 211 c c C REFLECTION ZONE SPECIFIED C PROBABLY THIS SHPULD BE PARSED -H +K +L = 3N : L=2N+1 : ETC C BUT IN THE SHORT TERM I WILL JUST SPECIFY 5 NUMBERS: C eg for above eqn -1 1 1 3 0 c for L=2N+1 0 0 1 2 1 2190 continue IRZ=IRZ+1 CALL GTPint(2,irzone(irz,1),NTOK,ITYP,FVALUE) CALL GTPint(3,irzone(irz,2),NTOK,ITYP,FVALUE) CALL GTPint(4,irzone(irz,3),NTOK,ITYP,FVALUE) CALL GTPint(5,irzone(irz,4),NTOK,ITYP,FVALUE) CALL GTPint(6,irzone(irz,5),NTOK,ITYP,FVALUE) go to 211 c Cc C NRFACTOR - Number of R factors to be printed... defaults to 15 2200 continue CALL GTPint(2,nrfact,NTOK,ITYP,FVALUE) if(nrfact.gt.500)then write(6,*) ' *** Only 500 R factors can be listed ****' nrfact=500 endif go to 211 c C LABI C C---- LABIN (File_Number) "a=b c=d ...." C ========== C 2175 continue 125 CONTINUE ITOK = 2 CALL LKYIN(MTZIN,LSPRGI,NLPRGI,NTOK,LINE,IBEG,IEND) go to 211 C 2170 continue Icorr=1 c if(nsymp.lt.2) call ccperr(1, 1 ' *** No translation search in P1 cell!!') IF(FMAX.EQ.0) FMAX=10000000. C BUILD UP SIN COS TABLE C PIBY2=6.28318 CONV=PIBY2/360. CONV1=8192./360. STEP=PIBY2/8192. THETA=0. DO 1110 I=1,2049 SINA(I)=SIN(THETA) THETA=THETA+STEP SINA(I+4096)=-SINA(I) SINA(8192+I)=SINA(I) SINA(4098-I)=SINA(I) SINA(8194-I)=-SINA(I) 1110 CONTINUE WRITE(6,3214) 3214 FORMAT(' SINS DONE') C DO 1101 I=1,8193,1024 C WRITE(6,3215)I,SINA(I),COSA(I) C1101 CONTINUE C3215 FORMAT(' CHECK SINE COS ',I5,2F10.5) C NFCT=0 SGFOOT=0. SIGFOT=0. DO 1100 I=1,5 NFC(I)=0 1100 SIGFO(I)=0.0000000 DO 1101 I=1,IGRDLM DELFT(I)=0. DO 1101 II=1,5 sigfco(I,ii)=0. sigfc(I,II)=0. sigfcc(I,II)=0. DELF(I,II)=0. 1101 CONTINUE C C C C open mtz file - extract cell dimensions first c cejd NLOOK=12+10*NSYM C C IPRINT = 2 C ======================================= CALL LROPEN(MTZIN,'HKLIN',IPRINT,IFAIL) C ======================================= CALL LRASSN(MTZIN,LSPRGI,NLPRGI,LOOKUP,CTPRGI) C ====================================================== CALL LRCELL(1,CELLMTZ) CALL LRNCOL(1,NCOLS) C IF (IFAIL.NE.0) CALL CCPERR(1, &'*** ERROR OPENING INPUT MTZ FILE ***') C CALL LRSYMI(1,NSYMP,LTYPE,NSPGRP,NAMSPG,PGNAM1) CALL LRSYMM(1,NSYM,RSYM1) C Chck that CAD symmetry had no tranlational elements, and C that the order of symops is the same. ILEN1 = LENSTR(PGNAM1) ILEN2 = LENSTR(PGNAM2) IF(PGNAM1(1:ILEN1).EQ.PGNAM2(1:ILEN2)) GO TO 1210 WRITE(6,*)' Error in symmetry - CADMTZ outputs PGNAME', 1 PGNAM1,' RSEARCH REQUESTS PGNAME',PGNAM2 1210 CONTINUE DO 1201 ISYM = 1,NSYMP DO 1202 I= 1,3 C This should no longer be necessary - I hope CAD does not apply C translational symmetry to partial phases. C IF(ABS(RSYM1(I,4,ISYM)).LE.0.0001 ) GO to 1211 C WRITE(6,*)' Error in symmetry - CADMTZ needs a', C 1 ' purely rotational spacegroup - eg PG2 PG3 PG422' C1211 CONTINUE DO 1203 J=1,3 IF( ABS(RSYM2(I,1,ISYM) - RSYM1(I,1,ISYM)) .LE.0.0001) 1 GO TO 1203 WRITE(6,*)' Error in symmetry - CAD and RSEARCH ', 1 'need symmetry operators in the same order - check SYMOP.LIB' call ccperr(1, 'Symmetry error') 1203 CONTINUE 1202 CONTINUE 1201 CONTINUE C ================================== C CALL SRLCF1(11,'HKLIN',-NLOOK,LOOK,LOOKUP,.TRUE.,NCOLS,CELLMTZ) C---- Now I can work out centricity and epsilon classes C C Centric refln defn. C ================================== IPRINT = 0 CALL CENTRIC(NSYM,RSYM1,IPRINT) CALL EPSLN (NSYM,NSYMP,RSYM1,IPRINT) C get ready for phase shifts. CAD output will not have done any. C Use input symmetry - not the rotational symmetry from CAD. C ********************************************************* CALL ASUSET(NAMSPG,NSPGRP,PGNAM2,NSYM,RSYM2,NSYMP,NLAUE,.TRUE.) C ********************************************************* C WRITE(6,1001)TITLE 1000 FORMAT(80A1) 1001 FORMAT(1X,A80) C C c WRITE(6,103) SCALE,BFACT,SMIN,SMAX,NSYMP,NSPGRP, 1 (NS,((rsym2(I,J,NS),J=1,4),I=1,4),NS=1,NSYMP) WRITE(6,203)FMIN,FMAX,FMAXC,map 103 FORMAT( 1 /' FP(USED) IS ',F10.5,'* FP(INPUT)*EXP(-', 1 F8.3,'(SIN THETA/LANDA)**2', 1 //,' 4S**2 LIMITS ARE ',2F10.5, 1 //,' NO OF SYMMETRY OPERATIONS ARE ',I5, ' spacegroup',i5, 1 48(//,' SYM OP ',I3,4X,4F8.3,3(/,15X,4F8.3))) 203 FORMAT( 1 //' REFLECTIONS EXCLUDED FROM R FACTOR SEARCH IF FOBS LT ',F10.5, 1 ' OR GT ',F15.2,//,' OR IF ANY FC TERM LT ',F10.5,//, 1 ' map section axis- ',i4,/, 1 ' ** If sec gt 0 "r factor" or "correlation" maps output ',/, 1 ' Section axis can be set but not fast and medium axes!!') if(ibin.gt.0.and.icorr.eq.0)write(6,*) 1 ' Rfactor map written to mapfile "mapout" ' if(ibcd.gt.0 .and. icorr.eq.0)write(6,*) 1 ' Rfactor maps listed to lineprinter' c if(ibin.gt.0.and.icorr.ne.0)write(6,*) 1 ' Fo Fc correlation map written to mapfile "mapout" ' if(ibcd.gt.0 .and. icorr.ne.0)write(6,*) 1 ' Fo Fc correlation map listed to lineprinter' c write(6,1113)nrfact 1113 format(//,' *** Number of R factors listed ***',i5) IF(BGRVAL.GT.0.01) 1 WRITE(6,1103)BGRVAL 1103 FORMAT(//,' BACKGROUND R VALUE SET TO ',F10.3) IF(BGRVAL.LT.0.01)WRITE(6,1104) 1104 FORMAT(//,' BACKGROUND R VALUE TAKEN AS AVERAGE R VALUE ') C C IF(ICEN.GT.0) WRITE(6,123) 123 FORMAT(///,' ONLY CENTRIC ZONES USED ') C IF(irz.GT.0) WRITE(6,124)Irz,((irzone(ir,j),j=1,5),ir=1,irz) 124 FORMAT(///,' Number of special refection zones',i5, 1 10(/,' Zone : ',i2,'*H + ',i2,'*K +',i2,'*L = ',i2,'n +',i2)) c C OPEN SEARCHSAVE - STATUS NEW IF IRESTT = 0 IF(IRESTT.EQ.0)CALL CCPOPN(22,'SEARCHSAVE',4,2,LDUM,IFAIL) c IF(IRESTT.GT.0)THEN C OPEN SEARCHSAVE - STATUS OLD CALL CCPOPN(22,'SEARCHSAVE',3,2,LDUM,IFAIL) WRITE(6,1105)NREF,IhH,IKk,ILl 1105 FORMAT(///,' RESTART AFTER ',I4,' REFLECTIONS.',/, 1 ' FIRST H K L IS ',3I4) c READ(22,END=1107)NREF1,IN,SIGFOT,SIGFO,sigfoo,sigfcc,sigfc,sigfco, 1 DELFT,DELF,NFCT,NFC IF(NREF1.eq.NREF) GO TO 1108 1107 call ccperr(1, 1 ' **** Something wrong with your Restart parameters**') c 1108 CONTINUE write(6,1109)nref1,in,sigfot 1109 format(' Check restart nref1 in sigfot..',i5,3x,3i4,f15.4) endif C extract cell dimensions to put in 'old' arrays AX=1./NX CELLM(1)=CELLMTZ(1) CELLM(4)=CELLMTZ(4) AAX=AX*CELLMTZ(1) WRITE(6,104)CELLMTZ(1),AAX,NX,IXMIN,IXMAX 104 FORMAT(/,' X CELL DIMENSION ',F6.2, 1 ' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - XLIMITS -',I5,' TO',I5) CELLM(2)=CELLMTZ(2) CELLM(5)=CELLMTZ(5) AY=1./NY AAY=AY*CELLMTZ(2) CELLM(3)=CELLMTZ(3) CELLM(6)=CELLMTZ(6) AZ=1./NZ AAZ=AZ*CELLMTZ(3) IF(ICODEY.EQ.0)WRITE(6,105)CELLMTZ(2),AAY,NY,IYMIN,IYMAX IF(ICODEY.EQ.1)WRITE(6,107)CELLMTZ(2),AAY,NY,IYMIN IF(ICODEY.EQ.2)WRITE(6,108)CELLMTZ(2),AAY,NY,IYMAX IF(ICODEY.EQ.3)WRITE(6,109)CELLMTZ(2),AAY,NY IF(ICODEZ.EQ.0)WRITE(6,106)CELLMTZ(3),AAZ,NZ,IZMIN,IZMAX IF(ICODEZ.EQ.4)WRITE(6,110)CELLMTZ(3),AAZ,NZ,IZMIN IF(ICODEZ.EQ.5)WRITE(6,111)CELLMTZ(3),AAZ,NZ,IZMIN IF(ICODEZ.EQ.6)WRITE(6,112)CELLMTZ(3),AAZ,NZ,IZMAX IF(ICODEZ.EQ.7)WRITE(6,113)CELLMTZ(3),AAZ,NZ,IZMAX IF(ICODEZ.EQ.8) WRITE(6,114)CELLMTZ(3),AAZ,NZ IF(ICODEZ.EQ.9)WRITE(6,115)CELLMTZ(3),AAZ,NZ 105 FORMAT(/ ' Y CELL DIMENSION ',F6.2, 1 ' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - YLIMITS -',I5,' TO',I5) 106 FORMAT(/ ' Z CELL DIMENSION ',F6.2, 1 ' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - ZLIMITS -',I5,' TO',I5) 107 FORMAT(/ ' Y CELL DIMENSION ',F6.2, 1' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - YLIMITS -',I5,' TO X' 1) 108 FORMAT(/ ' Y CELL DIMENSION ',F6.2, 1' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - YLIMITS -X TO',I5) 109 FORMAT(/ ' Y CELL DIMENSION ',F6.2, 1' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - YLIMITS -X TO X') 110 FORMAT(/ ' Z CELL DIMENSION ',F6.2, 1' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - ZLIMITS -',I5,' TO X' 1) 111 FORMAT(/ ' Z CELL DIMENSION ',F6.2, 1' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - ZLIMITS -',I5,' TO Y' 1) 112 FORMAT(/ ' Z CELL DIMENSION ',F6.2, 1' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - ZLIMITS -X TO',I5) 113 FORMAT(/ ' Z CELL DIMENSION ',F6.2, 1' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - ZLIMITS -Y TO',I5) 114 FORMAT(/ ' Z CELL DIMENSION ',F6.2, 1' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - ZLIMITS -X TO X') 115 FORMAT(/ ' Z CELL DIMENSION ',F6.2, 1' STEP SIZE(IN AS) ',F7.3,'(=1/',I3,') - ZLIMITS -Y TO Y') CC C c ejd 3-12-87 C BILD AN ARRAY IGRID(..) =2**20*IXPTS+2**10*IYPTS+IZPTS c FOR EACH GRID POINT . C C C ixpts=ixmax-ixmin+1 ICOUNT=1 npts=0 DO 8 IX1=1,IXpts IGR=ITWO20*(IX1-1) IF(ICODEY.EQ.1.OR.ICODEY.EQ.3) then iymin=ix1+ixmin-1 IF(ICODEY.EQ.2 .OR. ICODEY.EQ.3) THEN iymax=ix1+ixmin-1 endif endif iypts=iymax-iymin+1 CC c Z limits now IF(ICODEZ.EQ.4.or.icodez.eq.8) then izmin=ix1+ixmin-1 IF(ICODEZ.EQ.5 .OR. ICODEZ.EQ.8) THEN izmax=ix1+ixmin-1 endif endif c DO 8 IY1=1,iypts c Z limits now c IF(ICODEZ.EQ.6.or.icodez.eq.9) then izmin=iy1+iymin-1 IF(ICODEZ.EQ.7 .OR. ICODEZ.EQ.9) THEN izmax=iy1+iymin-1 endif endif c izpts=izmax-izmin+1 npts=npts + izpts if(npts.gt.igrdlm) then WRITE(6,1005) NPTS, IGRDLM 1005 FORMAT(//,' ', I8, ' grid points requested and limit is ',I8) call ccperr(1, + 'Set logical name RSEARCH_NSIZE to the size required') endif c IGRI=IGR+ITWO10*(IY1-1) DO 9 IZ1=1,izpts IGRID(ICOUNT)=IGRI+IZ1-1 ICOUNT=ICOUNT+1 9 CONTINUE c 8 CONTINUE c IF(IXpts.GT.ITWO10.OR.IYpts.GT.ITWO10.OR.IZpts.GT.ITWO10) * call ccperr(1, 'Need Grid Limits < 1024 - think again') C ICNTOT=ICOUNT-1 if(nrfact.gt.icntot)nrfact=icntot WRITE(6,1006)ICNTOT,cellm 1006 FORMAT(//,' Number of grid points to inspect is ',I10,/, 1 ' cell dimensions ',6f7.3) CC READ P1 SFS FILE FOR DIFFERENT SYMMETRY ORIENTATIONS CC C C READ FILE NAME AND DETAILS C 7 continue C C START TO READ H K L ETC C 510 CONTINUE C CALL LRREFF(MTZIN,SSQ,ADATA,EOF) IF(EOF) GO TO 530 CALL LRREFM(MTZIN,LOGMSS) C Require Fobs and SIGFobs to exist IF( LOGMSS(4) .OR. LOGMSS(5) ) GO TO 510 IN(1)=NINT(ADATA(1)) IN(2)=NINT(ADATA(2)) IN(3)=NINT(ADATA(3)) IF(IRESTT.GT.0)THEN IF(IHh.NE.IN(1)) GO TO 510 IF(IKk.NE.IN(2)) GO TO 510 IF(ILl.NE.IN(3)) GO TO 510 IRESTT=0 nref=nref-1 ENDIF c if (irz.gt.0) then do 511 ir=1,irz ichk= irzone(ir,1)*iN(1)+irzone(ir,2)*iN(2) 1 +irzone(ir,3)*iN(3) If(irzone(ir,4) .ne. 0) 1 ichk=mod(ichk+1000*irzone(ir,4),irzone(ir,4)) if(ichk.ne.irzone(ir,5) ) go to 510 511 continue endif c FO=0 SD=0 FC=0 C FO=ADATA(4) SD=ADATA(5) scb=EXP(-BFACT*SSQ/4.)*SCALE FO=FO*SCB IF(SD.LE.0.0000001) go to 510 IF(SSQ.LT.SMIN.OR.SSQ.GT.SMAX) go to 510 IF(ICEN.EQ.0) GO TO 210 C C CHECK FOR CENTRIC ZONES C call centr(in,ic) if(ic.ne.0) go to 210 go to 510 C 210 CONTINUE c Split the data into 5 resolution bins for statistics ASIN=5.*(SSQ-SMIN)/(SMAX-SMIN)+1.000001 ISIN=ASIN IF(ISIN.GT.5) ISIN=5 C IF(FO.GT.FMAX)THEN IXFPMX(ISIN)=IXFPMX(ISIN)+1 go to 510 ENDIF c IF(FO.LT.FMIN)THEN IXFPMN(ISIN)=IXFPMN(ISIN)+1 go to 510 ENDIF C c scb=EXP(-BFACT*SSQ/4.)*SCALE/10000.0 C C FORM PARTIAL A B TERMS FOR EACH SYMMETRY OPERATION c FCMIN=10000. c c loop round FCparts for each symmetry operation. c As long as ADATA(6) ( FC1 ) is there all others should be.. IF (LOGMSS(6)) GO TO 510 DO 11 ISYM=1,NSYMP IFC=6+2*(ISYM-1) IF (LOGMSS(IFC) .OR. LOGMSS(IFC+1)) THEN WRITE(6,*) 'This should not happen!! Missing FCx, PCx' GO TO 510 ENDIF PHASIN = ADATA(IFC+1) CALL ASUPHP(IN,ISYM,1,PHASIN,PHASOUT) ALPHA=CONV1*PHASOUT IALPH=ALPHA+1 IF(IALPH.LT.0) IALPH=IALPH+8192 FC=ADATA(IFC) c FC=scb*FC FC=FC IF(FC.LT.FCMIN)FCMIN=FC ACPART(ISYM)=FC*COSA(IALPH) BCPART(ISYM)=FC*SINA(IALPH) c c THIS Is Correct _ DONOT multiply by Rsymt C 6-2-89 OH bother - it wasn't... the indices of ROT are wrong . c It only affects 3 or 6 fold symmetry c c wrong IH =IN(1)*rsym(1,1,I)+IN(2)*rsym(1,2,I)+IN(3)*rsym(1,3,I) c wrong IK =IN(1)*rsym(2,1,I)+IN(2)*rsym(2,2,I)+IN(3)*rsym(2,3,I) c wrong IL =IN(1)*rsym(3,1,I)+IN(2)*rsym(3,2,I)+IN(3)*rsym(3,3,I) IH =IN(1)*rsym1(1,1,ISYM)+IN(2)*rsym1(2,1,ISYM) 1 +IN(3)*rsym1(3,1,ISYM) IK =IN(1)*rsym1(1,2,ISYM)+IN(2)*rsym1(2,2,ISYM) 1 +IN(3)*rsym1(3,2,ISYM) IL =IN(1)*rsym1(1,3,ISYM)+IN(2)*rsym1(2,3,ISYM) 1 +IN(3)*rsym1(3,3,ISYM) C Add phase shift..... c C NO NEED TO MULTIPLY BT PIBY2 WHEN USING LOOKUP TABLE FOR SIN/COS c Store ah ak al for each symmetry operation c AH(ISYM)=AX*IH AK(ISYM)=AY*IK AL(ISYM)=AZ*IL C 11 CONTINUE 13 CONTINUE C IF(FCMIN.LT.FMAXC)THEN IXFCMN(ISIN)=IXFCMN(ISIN)+1 go to 510 ENDIF C C ALL REJECTION TESTS PASSED - COUNT REFLECTIONS AND SAVE SUMS EACH 100 C REFLECTIONS C SAVE DELFT DELF SIGFOT SIGFO H K L NREF ETC FOR RESTART NREF=NREF+1 ICHK=MOD(NREF,100) IF(ICHK.EQ.1)THEN LPRINT = .TRUE. DO 133 jj = 4,NCOLS IF(LOGMSS(JJ)) LPRINT=.FALSE. 133 CONTINUE IF(LPRINT) WRITE(6,3459)NREF,IN,ssq,(adata(jj),jj=4,ncols) 3459 FORMAT(' REFLECTION ',I5,/,3X,3I4,f6.3,10f8.1,2(/,21x,10f8.1)) rewind (22) WRITE(22)NREF,IN,SIGFOT,SIGFO,sigfoo,sigfcc,sigfc,sigfco, 1 DELFT,DELF,NFCT,NFC ENDIF C NFCT=NFCT+1 SGFOOT=SGFOOT+FO*FO SIGFOO(ISIN)=SIGFOO(ISIN)+FO*FO SIGFOT=SIGFOT+FO SIGFO(ISIN)=SIGFO(ISIN)+FO NFC(ISIN)=NFC(ISIN)+1 CC 3456 FORMAT(3I5) c Check - Do you want to add an FC vector onto your SEARCH vector?? c If so - set column FC=... then idata(31) will be gt 0 C If FC has been set initial values equal fc*cos(PHIC) FC*sin(phic) c AC1=0. BC1=0. IF ((.NOT. LOGMSS(30)) .OR. (.NOT. LOGMSS(31))) THEN IF (ABS(adata(30)).gt.0.001) THEN ALPHA=CONV1*ADATA(31) IALPH=ALPHA+1 IF(IALPH.LT.0) IALPH=IALPH+8192 FC=ADATA(30) c FC=scb*FC AC1=FC*COSA(IALPH) BC1=FC*SINA(IALPH) ENDIF ENDIF C c Unpack each grid point requested for trial translation... DO 12 ICOUNT=1,ICNTOT IXa=IGRID(ICOUNT)/ITWO20 IGRI=IGRID(ICOUNT)-IXa*ITWO20 ix= ixa + ixmin IF(ICODEY.EQ.1.OR.ICODEY.EQ.3) then iymin=ix endif IYA=IGRI/ITWO10 iy= iya + iymin c Z limits now IF(ICODEZ.EQ.4.or.icodez.eq.8) then izmin=ix endif IF(ICODEZ.EQ.6.or.icodez.eq.9) then izmin=iy endif IZ=IGRI-ITWO10*IYA + izmin c c AC=AC1 BC=BC1 DO 17 I=1,NSYMP C COS ARRAY GOES FROM 0 TO 8192 TO COVER 360 DEGREES. DEL=8192.*(AH(I)*IX+AK(I)*IY+AL(I)*IZ+512.) IDEL=DEL IDELM=MOD(IDEL,8192) +1 COSDEL=COSA(IDELM) SINDEL=SINA(IDELM) AC=ACPART(I)*COSDEL-BCPART(I)*SINDEL + AC BC=ACPART(I)*SINDEL+BCPART(I)*COSDEL + BC 17 CONTINUE FCC=AC*AC+BC*BC FC=0.0 if(FCC.gt.0.0)Fc=sqrt(FCC) DELF(ICOUNT,ISIN)=DELF(ICOUNT,ISIN)+ABS(FO-FC) sigfcc(ICOUNT,ISIN)=sigfcc(ICOUNT,ISIN)+FCC SIGFC(ICOUNT,isin)=SIGFC(ICOUNT,isin)+FC sigfco(ICOUNT,isin)=sigfco(ICOUNT,isin)+FC*FO DELFT(ICOUNT)=DELFT(ICOUNT)+ABS(FO-FC) 12 CONTINUE go to 510 CC CC TABULATE RESULTS CC 530 WRITE(6,3496) 3496 FORMAT(' HKL INPUT DONE ') if(sigfot.eq.0.0) 1 call ccperr(1, '**** No reflections selected ***') DELTOT=0. ATOT=0. DO 331 I=1,5 331 DELTSN(I)=0. DO 31 ICOUNT=1,ICNTOT DELFT(ICOUNT)=100.*DELFT(ICOUNT)/SIGFOT DELTOT=DELTOT+DELFT(ICOUNT) DO 31 ISIN=1,5 if(sigfo(isin).gt.0.0) 1 DELF(ICOUNT,ISIN)=100.*DELF(ICOUNT,ISIN)/SIGFO(ISIN) DELTSN(ISIN)=DELTSN(ISIN)+DELF(ICOUNT,ISIN) 31 CONTINUE ATOT=ICNTOT DO 332 I=1,5 DELTSN(I)=DELTSN(I)/ATOT 332 CONTINUE DELTOT=DELTOT/ATOT WRITE(6,1007)ICNTOT,DELTOT,(DELTSN(I),I=1,5), 1 (IXFPMN(I),IXFPMX(I),IXFCMN(I),I=1,5) 1007 FORMAT(' AVERAGE R VALUE FOR ',I5,' GRID POINTS ',F10.5,/, 1 ' AVERAGE FOR EACH 4(S/L)**2 RANGE IS ',5(F10.5,5X),/, 1 ' EXCLUDED DATA - < FPMIN - > FPMAX - FC< FCMIN',5(/,12x,3I10)) DO 7777 I=1,nrfact FMIN=10000. JMIN=0 DO 7776 J=1,ICNTOT IF(FMIN.LE.DELFT(J)) GO TO 7776 JMIN=J FMIN=DELFT(J) 7776 CONTINUE FKP(I)=DELFT(JMIN) JKP(I)=JMIN DELFT(JMIN)=10000. 7777 CONTINUE CC CC WRITE(6,1008) 1008 FORMAT(//,' OUTPUT FOR LOWEST R VALUE POSITIONS IS',/, 1 ' 4(S/L)**2 LIMITS R% NREFLNS MEAN' 1 ,2X,'RMS SCALE ',//, 1 ' BEWARE !! Translations are ALONG CRYSTAL AXES -' , 1 ' NOT ORTHOGONAL AXES!!!!!') DO 37 I=1,nrfact ICOUNT=JKP(I) IXa=IGRID(ICOUNT)/ITWO20 IGRI=IGRID(ICOUNT)-IXa*ITWO20 ix= ixa + ixmin IF(ICODEY.EQ.1.OR.ICODEY.EQ.3) then iymin=ix endif IYA=IGRI/ITWO10 iy= iya + iymin c Z limits now IF(ICODEZ.EQ.4.or.icodez.eq.8) then izmin=ix endif IF(ICODEZ.EQ.6.or.icodez.eq.9) then izmin=iy endif IZ=IGRI-ITWO10*IYA + izmin C WRITE(6,*)IGRID(ICOUNT),IXA,IYA,IGRI,IX,IY,IZ,IXMIN,IYMIN,IZMIN c c IX=IGRID(ICOUNT)/ITWO20 c IGRI=IGRID(ICOUNT)-ITWO20*IX c IY=IGRI/ITWO10 c IZ=IGRI-ITWO10*IY c XSHIFT=AX*IX YSHIFT = AY*IY ZSHIFT=AZ*IZ XSHIFA=CELLMTZ(1)*AX*IX YSHIFA=CELLMTZ(2)* AY*IY ZSHIFA=CELLMTZ(3)*AZ*IZ DELFT(ICOUNT)=FKP(I) SFCFCO=sigfco(ICOUNT,1)+sigfco(ICOUNT,2)+sigfco(ICOUNT,3) 1 +sigfco(ICOUNT,4)+sigfco(ICOUNT,5) SFCFCT=sigfcc(ICOUNT,1)+sigfcc(ICOUNT,2)+sigfcc(ICOUNT,3) 1 +sigfcc(ICOUNT,4)+sigfcc(ICOUNT,5) FOMEAN=SIGFOT/NFCT rmsfc=0.0 if(sfcfct.gt.0.0)RMSFC=SQRT(SFCFCT/NFCT) WRITE(6,1002) 1 XSHIFT,YSHIFT,ZSHIFT,XSHIFA,YSHIFA,ZSHIFA,IX,IY,IZ 1002 FORMAT(///,' R FACTORS FOR TRANSLN ',3F7.4,'(FRACT)',/, 1 24X,3F7.2,'(IN AS ALONG CRYSTAL AXES)', 1 /,24X,3I7,'( GRID POINT)') SCAL=sfcfco/SFCFCT WRITE(6,1003)SMIN,SMAX,FKP(I),NFCT,FOMEAN,RMSFC,SCAL C C SMN=SMIN SMX=SMIN STEP=(SMAX-SMIN)/5. DO 38 ISIN=1,5 SMN=SMX SMX=SMX+STEP IF(NFC(ISIN).EQ.0) GO TO 38 FOMEAN=SIGFO(ISIN)/NFC(ISIN) rmsfc=0.0 if(sigfcc(icount,isin).gt.0.0) 1 RMSFC=SQRT(sigfcc(ICOUNT,ISIN)/NFC(ISIN)) WRITE(6,1003)SMN,SMX,DELF(ICOUNT,ISIN),NFC(ISIN),FOMEAN,RMSFC 38 CONTINUE 1003 FORMAT(2F8.4,5X,F10.2,5X,I5,5X,2F10.2,F10.5) 37 CONTINUE C C PRINT R FACTOR MAP IF REQUIRED - section axis chosen to get fewest c sections. No map unless grid is "rectangular", ie icodey and icodez =0 c if(icodey.gt.0 .or.icodez.gt.0) then write(6,*)' ** No map output - irregular shape***' map=0 endif IF (MAP.EQ.0) GO TO 39 IDIG=4 ISEC=1 ILIN=1 ISP=0 IPERM=0 VOL=1 F000=0 itype=-1 NCOL=29 NLINES=60 PERM(1,1)=1 PERM(2,2)=1 PERM(3,3)=1 PERM(4,4)=1 c c set kp1 2 3 for map requirements c IF(MAP.EQ.2)then nf=izpts nm=ixpts nsec=iypts nsec0=iymin kp1=1 kp2=2 kp3=3 Iscpts=ixpts*izpts endif c IF(MAP.EQ.3)then nf=ixpts nm=iypts nsec=izpts nsec0=izmin kp1=2 kp2=3 kp3=1 Iscpts=ixpts*iypts endif c IF(MAP.EQ.1)then nf=iypts nm=izpts nsec=ixpts nsec0=ixmin kp1=3 kp2=1 kp3=2 Iscpts=iypts*izpts endif c if(iscpts.gt.jden)WRITE(6,1011)JDEN,Iscpts 1011 FORMAT(///,' MAP SECTION TOO LARGE - LIMIT IS ',I5, 1 ' NEED ',I6) IF(Iscpts.GT.JDEN) GO TO 39 c c SMX=SMIN c c Six maps output to lineprinter if requested c - one for all data and one each for each resln range. c c Only the complete one output to "mapout" c DO 441 IMAP=1,6 C WE CAN INPUT A BACKGROUND R VALUE C c this is to calculate the CORRELATION COEFFICIENT of two variables c ****************************************************************** c ** ave(X*Y) - ave(X)*ave(Y) c * corr. coeff.= ------------------------------------------------- c ** SQRT(ave(X**2)-ave(X)**2)*SQRT(ave(Y**2)-ave(Y)**2) c ****************************************************************** BGR=BGRVAL BGRS=BGRVAL IF(IMAP.EQ.1.AND.BGR.LT.0.01)BGR=DELTOT IF(IMAP.GT.1.AND.BGRS.LT.0.01)BGRS=DELTSN(ISIN-1) c c map for all data Iyzpts=iypts*izpts IF(IMAP.EQ.1) GO TO 44 c ibin=0 if(ibcd.eq.0) go to 39 ISIN=IMAP-1 SMN=SMX SMX=SMN+STEP IF(NFC(ISIN).EQ.0) GO TO 441 c 44 continue c Calculate avfoo avfo once only IF(icorr.gt.0)then c IF(IMAP.EQ.1 )then avfoo= (sigfoo(1) + sigfoo(2) + sigfoo(3) + sigfoo(4) + sigfoo(5)) 1 /nfct avfo= (sigfo(1) + sigfo(2) + sigfo(3) + sigfo(4) + sigfo(5)) 1 /nfct endif c IF(IMAP.gt.1 )then avfoo= sigfoo(isin)/nfc(isin) avfo= sigfo(isin)/nfc(isin) endif ajunk1 = (avfoo - avfo*avfo) if(ajunk1.le.0.0) 1 call ccperr(1, 1 ' Something very wrong - (avfoo - avfo*avfo) lt 0') ajunk1=sqrt(ajunk1) endif DO 41 Isec=nsec0,nsec0 + nsec-1 DO 42 IXY=1,iscpts 42 DEN(IXY)=0. c write(6,*) ' MAP ISEC nm nf iscpts iyzpts,izpts ', c 1 MAP,ISEC,nm,nf,iscpts,iyzpts,izpts C JD=0 C DO 43 Im=1,nm DO 43 Iff=1,nf C MAP=2 ISEC Y NF Z NM X IF(MAP.EQ.2)ITEM = (IM-1)*IYZPTS + (ISEC-nsec0)*IZPTS + IFF C MAP=3 ISEC Z NF X NM Y IF(MAP.EQ.3) c ITEM = (IFF-1)*IYZPTS + (IM-1)*IZPTS + ISEC-nsec0 + 1 C MAP=1 ISEC X NF Y NM Z IF(MAP.EQ.1)ITEM = (ISEC-nsec0)*IYZPTS + (IFF-1)*IZPTS + IM c JD = JD + 1 c this is to calculate the CORRELATION COEFFICIENT of two variables c ****************************************************************** c ** ave(X*Y) - ave(X)*ave(Y) c * corr. coeff.= ------------------------------------------------- c ** SQRT(ave(X**2)-ave(X)**2)*SQRT(ave(Y**2)-ave(Y)**2) c ****************************************************************** FOMEAN=SIGFOT/NFCT IF(icorr.gt.0)then c IF(IMAP.EQ.1 )then avfco= (sigfco(item,1) + sigfco(item,2) + sigfco(item,3) + 1 sigfco(item,4) + sigfco(item,5) )/nfct avfc= (sigfc(item,1) + sigfc(item,2) + sigfc(item,3) + 1 sigfc(item,4) + sigfc(item,5) )/nfct avfcc= (sigfcc(item,1) + sigfcc(item,2) + sigfcc(item,3) + 1 sigfcc(item,4) + sigfcc(item,5) )/nfct endif c IF(IMAP.gt.1 )then avfco= sigfco(item,isin)/nfc(isin) avfcc= sigfcc(item,isin)/nfc(isin) avfc= sigfc(item,isin)/nfc(isin) avfcc= sigfcc(item,isin)/nfc(isin) endif c acor1= avfco - avfo*avfc c ajunk2 = (avfcc - avfc*avfc) if(ajunk2.le.0.0) 1 call ccperr(1, 1 ' Something very wrong - (avfcc - avfc*avfc) lt 0') ajunk2=sqrt(ajunk2) c acor2= ajunk1*ajunk2 den(jd) = acor1/acor2 endif if(icorr.eq.0)then DEN(JD)=BGR-DELFT(ITEM) IF(IMAP.GT.1)DEN(JD)=BGRS-DELF(ITEM,ISIN) endif 43 CONTINUE IF(IMAP.EQ.1.and.icorr.eq.0) WRITE(6,1009) 1009 FORMAT('2 MAP OF 10.*(BACKGROUND R VALUE - GRID R VALUE)', 1 ' - FULL SIN THETA RANGE ') IF(IMAP.GT.1.and.icorr.eq.0) WRITE(6,1010)SMN,SMX 1010 FORMAT('2 MAP OF 10.*(BACKGROUND R VALUE - GRID R VALUE)', 1 ' 4(S/L)**2 RANGE IS ',2F8.4) IF(IMAP.EQ.1.and.icorr.gt.0) WRITE(6,1019) 1019 FORMAT('2 MAP OF Correlation between Fo and Fc', 1 ' - FULL SIN THETA RANGE ') IF(IMAP.GT.1.and.icorr.gt.0) WRITE(6,1020)SMN,SMX 1020 FORMAT('2 MAP OF Correlation between Fo and Fc', 1 ' 4(S/L)**2 RANGE IS ',2F8.4) nsc=1 CALL PRNTYZ(DEN,NF,NM,NSC,Isec,IBIN,IBCD) 41 CONTINUE 441 CONTINUE 39 RETURN 99 WRITE(6,1004) 1004 FORMAT(' SOMETHING WRONG WITH LIMIT CARDS') RETURN END C SUBROUTINE PRNTYZ(X, Nf, Nm, Ns, isec0, BIN, BCD) C ============================================= C INTEGER nf,nm,ns,isec0, BIN, BCD c REAL X(NY,NX,NZ) REAL X(Nf,nm,ns) C C FFT map output subroutine for section defined by 'kp3' C This version writes the binary file with 'kp2' running fastest C LOGICAL SPEW C INTEGER XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX,xyzmn(3),xyzmx(3) INTEGER TYPE, NLINES, NCOL, IUVW(3) REAL V, F000, RHOMAX, RHOMIN C COMMON /ATSYM/rsym1(4,4,192),Rsym2(4,4,192),NSYM,PERM(4,4), 1 MATOMS,LATOMS COMMON /mphdr/ NXYZ(3),XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX, . NUSED6,NUSED7,NUSED9,NUSED10,NUSED11,IDUMMY(3),NUSED12 . ,CELL(6),LspGRP,RHMIN,RHMAX,Rhmean,rhrms INTEGER XX(6) EQUIVALENCE (XX,XMIN) COMMON /IOFFT/ T2,T3,T4,NSYMCH,NABCH,NSECAX COMMON /FFTSPG/NSPGRP CHARACTER*80 TITLE COMMON/mphdrr/TITLE COMMON /MAPFMT/ TYPE, NLINES,NCOL, V, F000, RHOMAX, RHOMIN COMMON/FOUT/MIN,MAX,ISP,IDIG,ILIN COMMON /FFTPRM/ KP1,KP2,KP3 INTEGER I, J, K, L, SEC, XL, XU, YZL, YZU, xllp,xulp,yzllp,yzulp INTEGER LINE(120) CHARACTER*1 KXYZ(3),AST,SPC REAL R, RMAX, RMIN CHARACTER*4 IFMT(8), JDIG(9),JLIN(4) SAVE DATA KXYZ/'X','Y','Z'/ DATA AST/'*'/,SPC/' '/ DATA IFMT/'(',' ','1X,I','3,1H','*,','128I','4',')'/ DATA JDIG/'1','2','3','4','5','6','7','8','9'/ DATA JLIN/' ','/ ','// ','/// '/ C c Take care.... nf nm ns are " nz nx p2" for FFT1 ,etc c AND nxyz(3) are the permuted values of NX NY NZ c SO nf and nm can be reset by PRMVECT and this can mess up the c definition of SPEW c nff=nf nmm=nm nss=ns C Now permute "nx" "ny" "nz" to get nxyz(..) related to A B C C c PRMVECT will return NXYZ AND XMIN.. to original order c CALL PRMVECT(PERM,NXYZ,1,1) CALL PRMVECT(PERM,XX(1),1,2) CALL PRMVECT(PERM,XX(2),1,2) C C IF THIS IS THE FIRST TIME, DERIVE THE SCALE FACTORS C xyzmn(1)=xmin xyzmn(2)=ymin xyzmn(3)=zmin C xyzmx(1)=xmax xyzmx(2)=ymax xyzmx(3)=zmax C C USE NSECAX TO CHOOSE WHETHER THIS IS A "PRINTY" LOOK-ALIKE C OR A "PRINTZ" ONE. C C FOR "PRINTZ" PUT kx1=kp1 kx2=kp2 kx3=kp3 kx1=kp1 kx2=kp2 kx3=kp3 C FOR "PRINTY" PUT kx1=kp1 kx2=kp3 kx3=kp2 if(nsecax.eq.2) then kx2=kp3 kx3=kp2 endif IF(isec0.NE.xyzmn(kx3)) GO TO 50 IF (NLINES .LE. 0) NLINES = xyzmx(kx1)-Xyzmn(kx1)+1 I=128/IDIG IF(NCOL.LE.0.OR.NCOL.GT.I)NCOL=I IF(BIN.LE.0) GO TO 1 C C Write out binary map file header C Title, NSEC (number of sections in map), fast, medium and sMIN axis C numbers, NNX,NNY,NNZ (sampling intervals along real xyz) C NSEC=xyzmx(kx3)-xyzmn(kx3)+1 C Set up map header in common C Axis order IUVW(1)=KX2 IUVW(2)=KX1 IUVW(3)=KX3 mode=2 CALL MWRHDR(BIN,TITLE,NSEC,IUVW,NXYZ,xyzmn(kx3), 1 xyzmn(kx2),xyzmx(kx2),xyzmn(kx1),xyzmx(kx1),CELL,LSPGRP,mode) C Copy symmetry operations from library file stream 4 to output map CALL MSYPUT(NSYMCH,LSPGRP,BIN) C c c XL XU refer to column numbers c YZL YZU refer to row numbers c 1 continue XL = xyzmn(kx1)+1 XU = xyzmx(kx1)+1 YZL = xyzmn(kx2)+1 YZU = xyzmx(kx2)+1 RMINOV=10000.0 RMAXOV=-10000.0 SPEW=(XL.EQ.1).AND.(XU.EQ.Nmm).AND. . (YZL.EQ.1).AND.(YZU.EQ.Nff).AND. (BIN.NE.0) IF(V.EQ.0.0) V = 1.0 IF (RHOMIN .LE. RHOMAX) GO TO 10 R = RHOMIN RHOMIN = RHOMAX RHOMAX = R 10 IF(TYPE.LT.0) GO TO 40 RMAX = X(YZL,XL,1) RMIN = RMAX c DO 20 I = 1,NZ DO 20 I = 1,Nss DO 20 J = XL,XU DO 20 K = YZL,YZU R = X(K,J,I) IF (R .GT. RMAX) RMAX = R IF (R .LT. RMIN) RMIN = R 20 CONTINUE IF(TYPE.GT.0) GO TO 30 RMAX = AMAX1(ABS(RMIN),RMAX) V = RMAX/RHOMAX F000 = 0.0 GO TO 40 30 V = (RMAX-RMIN)/(RHOMAX-RHOMIN) F000 = V*RHOMAX-RMAX 40 CONTINUE WRITE(6,1000) F000,V IF(TYPE.GE.0)WRITE(6,2000)Nss WRITE(6,6001) KXYZ(kx3) 6001 FORMAT(/' Section axis is ',A1) IF(BIN.NE.0) WRITE(6,6002) KXYZ(kx2),KXYZ(kx1) 6002 FORMAT(/' Binary output has ',A1,' running fastest, ',A1,' next') IF(BCD.NE.0) WRITE(6,6003) KXYZ(kx2),KXYZ(kx1) 6003 FORMAT(/' Line-printer output has ',A1,' across, ',A1, . ' down the page') C C Initialise max,min,mean RHMEAN=0.0 RHMIN=(X(YZL,XL,1)+F000)/V RHMAX=RHMIN RMS=0.0 NPT=0 C C PUT OUT SCALED "YZ" SECTIONS C 50 CONTINUE DO 700 K = 1,Nss SEC = isec0+K-1 C C SCALE THIS SECTION C XL = xyzmn(kx1)+1 XU = xyzmx(kx1)+1 YZL = xyzmn(kx2)+1 YZU = xyzmx(kx2)+1 RMIN = (F000+X(YZL,XL,K))/V RMAX = RMIN kx1mIN=XL-1 kx2mIN=YZL-1 kx1mAX=XL-1 kx2mAX=YZL-1 DO 100 J = XL,XU DO 100 I = YZL,YZU R = (F000+X(I,J,K))/V X(I,J,K) = R RHMEAN=RHMEAN+R RMS=RMS+R*R NPT=NPT+1 CCC IF(R.GT.RMAX)RMAX = R CCC IF(R.LT.RMIN)RMIN = R CCC100 CONTINUE C C FIND MIN AND MAX VALUES FOR THIS SECTION : FIND POSITION THEY OCCU C AND ALSO OVERALL MIN AND MAX. C IF(R.GE.RMIN) GO TO 301 RMIN=R kx1min=J-1 kx2min=I-1 301 IF(R.LE.RMAX) GO TO 100 RMAX=R kx1max=J-1 kx2max=I-1 100 CONTINUE c c Find overall max min IF(RMAX.LE.RMAXOV) GO TO 302 RMAXOV=RMAX kx1mxv=kx1max kx2mxv=kx2max kx3mxv=SEC 302 IF(RMIN.GE.RMINOV) GO TO 303 RMINOV=RMIN kx1mnv=kx1min kx2mnv=kx2min kx3mnv=SEC 303 CONTINUE RHMIN=AMIN1(RHMIN,RMIN) RHMAX=AMAX1(RHMAX,RMAX) RHMIN=AMIN1(RHMIN,RMIN) RHMAX=AMAX1(RHMAX,RMAX) C C BINARY OUTPUT C IF(BIN.EQ.0) GO TO 160 IF(SPEW) GO TO 150 cejd CALL MWRSEC(BIN,X(1,1,K),NY,NX,YZL,YZU,XL,XU) CALL MWRSEC(BIN,X(1,1,K),Nff,Nmm,YZL,YZU,XL,XU) GO TO 160 C 150 CALL MSPEW(BIN,X(1,1,K)) 160 CONTINUE C C DECIMAL OUTPUT C IF(BCD.EQ.0) GO TO 600 c Increase scale if rhomax lt 1 ichk =nint(100.*rhomax) if(ichk.lt.10)scc=100. if(ichk.lt.100)scc=10. if(ichk.ge.100)scc=1. if(scc.gt.1.)write(6,9753) scc 9753 format(' ***** Map scale multiplied by ',f6.2, 1 ' to get numbers > 1.0 *****') C C CONSTRUCT FORMAT IFMT(2)=JLIN(ILIN) IFMT(7)=JDIG(IDIG) YZUlp = xyzmn(kx2) 200 CONTINUE YZLlp = YZUlp+1 YZUlp = YZUlp + NCOL IF(YZUlp.GT.xyzmx(kx2)) YZUlp = xyzmx(kx2)+1 XUlp = xyzmn(kx1) 300 CONTINUE XLlp = XUlp+1 XUlp = XUlp+NLINES IF(XUlp.GT.xyzmx(kx1))XUlp = xyzmx(kx1)+1 I = YZLlp-1 J = YZUlp-1 WRITE(BCD,3000) TITLE,KXYZ(kx3),SEC,KXYZ(kx2),I,J L=0 DO550J=YZLlp,YZUlp L=L+1 550 LINE(L)=J-1 WRITE(BCD,IFMT)SEC,(LINE(J),J=1,L) N=IDIG-1 WRITE(BCD,4600)( (SPC,I=1,N), AST, J=1,L) DO 500 J = XLlp,XUlp LINE(1) = J-1 L = 1 DO 400 I = YZLlp,YZUlp L = L+1 LINE(L)=NINT(scc*X(I,J,K)) IF(LINE(L).GT.MIN.AND.LINE(L).LT.MAX)LINE(L)=0 400 CONTINUE WRITE(BCD,IFMT)(LINE(I),I = 1,L) 500 CONTINUE IF(XUlp.LE.xyzmx(kx1)) GO TO 300 IF(YZUlp.LE.xyzmx(kx2)) GO TO 200 C 600 CONTINUE CCC WRITE(6,5000)SEC,RMIN,RMAX IF(BCD.GT.0) WRITE(6,4999) WRITE(6,5000)KXYZ(kx3),SEC,RMIN,KXYZ(kx1),kx1min,KXYZ(kx2),kx2min, 1 RMAX,KXYZ(kx1),kx1max,KXYZ(kx2),kx2max 700 CONTINUE C IF ( ISEC0 + Nss .LE. xyzmx(kx3)) GOTO 800 C IF(BIN.eq.0) go to 880 C Close map file - MCLOSE prints out density limits C CALL MCLOSE(BIN,RHMIN,RHMAX,RHMEAN,RMS) 880 continue TMEAN2 = (RHMEAN/FLOAT(NPT))**2 RHRMS = SQRT(RMS/FLOAT(NPT) - TMEAN2) WRITE(6,5050) RHRMS C C Print max min information C WRITE(6,5001)RMAXOV,KXYZ(kx3),kx3mxv,KXYZ(kx1),kx1mxv, 1 KXYZ(kx2),kx2mxv RHMEAN = RHMEAN/FLOAT(NPT) WRITE(6,5002)RMINOV,KXYZ(kx3),kx3mnv,KXYZ(kx1),kx1mnv, 1 KXYZ(kx2),kx2mnv,RHMEAN C 800 CONTINUE c Restore permuted order of nxyz xmin etc CALL PRMVEC(PERM,NXYZ,1,1) CALL PRMVEC(PERM,XX(1),1,2) CALL PRMVEC(PERM,XX(2),1,2) RETURN C 1000 FORMAT('0F000 =',1PE12.3,10X,'V =',E12.3) 2000 FORMAT('0SCALING BASED ON VALUES IN THE FIRST',I3,' SECTIONS.') 3000 FORMAT(////10X,A80,5X,A1,' =',I3,', ',A1,' =',I3,' TO',I4/) 4600 FORMAT(5X,128A1) CC5000 FORMAT(10H0Section =,I4,19H minimum density =,F15.4, CC - 19H maximum density =,F15.4) 4999 FORMAT(1X) 5000 FORMAT(' Section ',A1,' = ',I3,' minimum rho =',F15.4, 1 ' at ',A1,' = ',I3,' ',A1,' = ',I3,' maximum rho =',F15.4, 2 ' at ',A1,' = ',I3,' ',A1,' = ',I3) 5001 FORMAT(//,' Overall maximum rho = ',F15.4, 1 ' on section ',A1,' = ',I3,' at ',A1,' = ',I3,1X,A1,' = ',I3) 5002 FORMAT(//,' Overall minimum rho = ',F15.4, 1 ' on section ',A1,' = ',I3,' at ',A1,' = ',I3,1X,A1,' = ',I3, 1 //,' Average density on this map is ',F19.8) 5050 FORMAT(30X,' Rms density =',F19.8) C END SUBROUTINE PRMVEC(PERM,JV,N,N1) C ============================= C C Permute vector JV(N,3) by permutation matrix PERM C N1 is first dimension of JV C DIMENSION JV(N1,3),LV(3),PERM(4,4) C C Permute DO 1 I=1,3 LV(I)=PERM(I,1)*JV(N,1)+PERM(I,2)*JV(N,2)+PERM(I,3)*JV(N,3) 1 CONTINUE C C Copy back DO 2 I=1,3 2 JV(N,I)=LV(I) RETURN END SUBROUTINE PRMVECT(PERM,JV,N,N1) C ============================= C C Permute vector JV(N,3) by transpose of permutation matrix PERM C N1 is first dimension of JV C DIMENSION JV(N1,3),LV(3),PERM(4,4) C C Permute DO 1 I=1,3 LV(I)=PERM(1,I)*JV(N,1)+PERM(2,I)*JV(N,2)+PERM(3,I)*JV(N,3) 1 CONTINUE C C Copy back DO 2 I=1,3 2 JV(N,I)=LV(I) RETURN END