PROGRAM FSFTOCTR C ================ C C ------------------------------------------------ C MAP REFORMATTING AND TRANSLATION PROGRAM LINKING C "FSFOR" STYLE MAPS TO "CONTUR" STYLE MAPS C BdV, 7-MAR-1989 C ------------------------------------------------ C C In the following example, we want a "CONTUR" map extending (in C fractional coordinates) from 0.333 to 1.333 in U, from 0.250 to C 1.500 in V, and from 0.750 to 1.250 in W (note that the CONTUR C map can cover more than one unit cell translation!) : C C input map output map C --------- ---------- C map direction: U V W U V W C grid : 48 64 80 48 64 80 C origin : 0 0 0 16 16 60 C last poin : 47 63 79 64 96 100 C #points : 48 64 80 49 83 41 C C C ** NOTE : the maximum number of points in one section of the C input FSFOR map is given in the following PARAMETER C statement. If this number must be changed, then change C all instances of this statement throughout the program! C C*** PARAMETER (MAPDIM=40000, MI=1, MO=2, LP=3, II=5, IO=6) C DIMENSION ITLE(20),XYZ(3),NAME(10),MAP(MAPDIM) COMMON /MAPS/ AMAP(MAPDIM) COMMON /GRID/ NXYZ(3), IORI(3), IEND(3), IOUT1(2), IOUT2(2) COMMON /ORTH/ YESNO EQUIVALENCE (MAP(1),AMAP(1)) CHARACTER*40, INFILE, OUTFILE, XYZ*1, ANS*1, NAME*10 LOGICAL YESNO DATA XYZ/'X','Y','Z'/, NTYP/1/ DATA NAME /'FSFOR',' ',' ',' ',' ', . ' ',' ',' ',' ',' '/ RMAX = -9999.9 RMS = 0.0 WRITE (IO,10) C C C-- MAP FILE NAMES C WRITE (IO,20) READ (II,90) INFILE OPEN (UNIT=MI,FILE=INFILE,STATUS='OLD',FORM='UNFORMATTED', . READONLY, ERR=9999) 1 WRITE (IO,30) READ (II,90) OUTFILE 2 CONTINUE MAPTYP = 1 CALL INITLP(LP,INFILE,OUTFILE,MAPTYP,NAME,NTYP) C C C-- READ IN HEADER TO GET INPUT MAP INFO C CALL HDRFSFOR(IO,LP,MI,IS,IA,ID) C C C-- REPORT INFO AND DETERMINE EXTENT OF OUTPUT MAP C 200 CONTINUE WRITE (IO,50) XYZ(IS),XYZ(IA),XYZ(ID),NXYZ(IS),NXYZ(IA),NXYZ(ID), . IORI(IS),IEND(IS),IORI(IA),IEND(IA),IORI(ID),IEND(ID) WRITE (LP,51) XYZ(IS),NXYZ(IS),IORI(IS),IEND(IS),XYZ(IA),NXYZ(IA), . IORI(IA),IEND(IA),XYZ(ID),NXYZ(ID),IORI(ID),IEND(ID) NPNTS = NXYZ(IA) * NXYZ(ID) IF (NPNTS.GT.MAPDIM) THEN WRITE (IO,60) NPNTS, MAPDIM WRITE (LP,61) NPNTS, MAPDIM STOP ENDIF CALL ORIENT(II,IO,IS,IA,ID,JS,JA,JD,NORIEN,IOU,IOV) OPEN (UNIT=MO,FILE=OUTFILE,STATUS='NEW',FORM='UNFORMATTED') C TYPE *, ' ' IOUT1(1) = IORI(JA) IOUT2(1) = IEND(JA) IOUT1(2) = IORI(JD) IOUT2(2) = IEND(JD) WRITE (6,75) XYZ(JA),IOUT1(1),IOUT2(1) READ (5,40) I,J IF (J.NE.0) THEN IOUT1(1) = I IOUT2(1) = J ENDIF WRITE (6,76) XYZ(JD),IOUT1(2),IOUT2(2) READ (5,40) I,J IF (J.NE.0) THEN IOUT1(2) = I IOUT2(2) = J ENDIF NF = IORI(JS) NL = IEND(JS) WRITE (6,77) XYZ(JS),NF,NL READ (5,40) I,J IF (J.NE.0) THEN NF = I NL = J ENDIF WRITE (IO,53) XYZ(IS),XYZ(IA),XYZ(ID),NXYZ(IS),NXYZ(IA),NXYZ(ID), . NF,NL,IOUT1(1),IOUT2(1),IOUT1(2),IOUT2(2) WRITE (LP,52) XYZ(JS),NXYZ(JS),NF,NL,XYZ(IA),NXYZ(IA), . IOUT1(1),IOUT2(1),XYZ(ID),NXYZ(ID),IOUT1(2),IOUT2(2) C C C-- POSITION MAP TO FIRST SECTION REQUIRED C CALL POSFSF(MI,0,NF,NXYZ(ID)) ICURR = NF C C C-- COPY THE REQUIRED SECTIONS TO THE OUTPUT FILE C DO ISEC=NF,NL IF (ICURR.GE.NXYZ(IS)) THEN ICURR = ICURR-NXYZ(IS) CALL POSFSF(MI,0,ICURR,NXYZ(ID)) ENDIF CALL RDMFSFOR(MI,IO,IS,IA,ID,RMAX,RMS) ICURR = ICURR+1 CALL WRMCTR(MO,IO,LP,OUTFILE,IS,IA,ID,JS,JA,JD,NF,NL, . NORIEN,RMAX,MAPTYP,ISEC) ENDDO C C C-- FINAL LINE ON THE CONTUR MAP C AOM7 = -1000.0 CCC = 999.0 / RMAX WRITE (MO) AOM7, RMAX, CCC C C C-- THAT'S IT: REPORT MAX AND RMS C CLOSE (UNIT=MO) NPNTS = NPNTS * (NL-NF+1) RMS = SQRT (RMS/NPNTS) WRITE (LP,81) RMAX, RMS STOP C C 9999 WRITE (IO,70) INFILE READ (II,90) INFILE OPEN (UNIT=MI,FILE=INFILE,STATUS='OLD',FORM='UNFORMATTED', . ERR=9999) GOTO 1 C 10 FORMAT (/T12,'F S F T O C T R --- V E R S I O N M A R C H', . ' 89'/) 20 FORMAT (T2,'Enter file name for the old (input) map',T54,'> ',$) 30 FORMAT (T2,'Enter file name for the new (output) map',T54,'> ',$) 40 FORMAT (2I5) 50 FORMAT(/T2,'The old map has ',A,' sections, ',A,' across and ', . A,' down'/T4,'with #grid :',T19,I4,T33,I4,T48,I4 . /T4,'from - to :',T17,I3,' -',I4,T31,I3,' -',I4,T46,I3,' -',I4) 53 FORMAT(/T2,'The new map has ',A,' sections, ',A,' across and ', . A,' down'/T4,'with #grid :',T19,I4,T33,I4,T48,I4,/T4, . 'from - to :',T17,I3,' -',I4,T31,I3,' -',I4,T46,I3,' -',I4/) 51 FORMAT (//T20,'DESCRIPTION OF INPUT FSFOR MAP'/ . /T58,'LAYOUT',T76,'GRID',T94,'FROM ... TO'//T40, . 'SECTIONS :',T61,A,T75,I4,T93,I4,5X,I4/T40,'ACROSS :',T61,A, . T75,I4,T93,I4,5X,I4/T40,'DOWN :',T61,A,T75,I4,T93,I4,5X,I4) 52 FORMAT (//T20,'DESCRIPTION OF OUTPUT CONTUR MAP'/ . /T58,'LAYOUT',T76,'GRID',T94,'FROM ... TO'//T40, . 'SECTIONS :',T61,A,T75,I4,T93,I4,5X,I4/T40,'ACROSS :',T61,A, . T75,I4,T93,I4,5X,I4/T40,'DOWN :',T61,A,T75,I4,T93,I4,5X,I4) 60 FORMAT (//T30,'*** MAP DIMENSION TOO LARGE ; CHANGE "MAPDIM" ***' . /T20,'(YOUR MAP HAS',I8,' POINTS, CURRENT DIMENSION IS',I8,')') 61 FORMAT (//T30,'*** MAP DIMENSION TOO LARGE ; CHANGE "MAPDIM" ***' . /T20,'(YOUR MAP HAS',I8,' POINTS, CURRENT DIMENSION IS',I8,')') 70 FORMAT (T2,'** Error opening file',T56,A . /T2,' Enter correct file name',T54,'> ',$) 75 FORMAT (' Enter first, last point in ',A,' (across)', . T42,'[',I4,',',I4,'] > ',$) 76 FORMAT (' Enter first, last point in ',A,' (down)', . T42,'[',I4,',',I4,'] > ',$) 77 FORMAT (' Enter first, last section (',A,')', . T42,'[',I4,',',I4,'] > ',$) 80 FORMAT (T2,'Maximum and rms density :',10X,F10.4,2X,F9.4/) 81 FORMAT (/T20,'MAXIMUM AND RMS DENSITY :',10X,F10.4,5X,F9.4) 90 FORMAT (A) END C C========================================================================== C SUBROUTINE INITLP(LP,INFILE,OUTFILE,MAPTYP,NAME,NTYP) DIMENSION NAME(NTYP) CHARACTER*10 DAT, TIM, NAME, INFILE*14, OUTFILE*14 WRITE (LP,10) INFILE, OUTFILE, MAPTYP, NAME(MAPTYP) RETURN C 10 FORMAT (///T44,'P R O G R A M F S F T O C T R'//T20, . 'READS ELECTRON DENSITY MAPS AND OUTPUTS THEM IN' . ' CONTUR FORMAT'//////T20, . 'THE INPUT MAP FILE IS ',A,5X,'THE OUTPUT MAP WILL BE WRITTEN ', . 'TO ',A//T20,'MAPTYPE IS ',I3,' : INPUT MAP FROM THE PROGRAM ' . ,A//T20,'DESCRIPTION OF THE OLD MAP :') END C C=========================================================================== C SUBROUTINE HDRFSFOR(IO,LP,MI,IS,IA,ID) C*** C This subroutine reads the header records on a FSFOR map file : C 1. a record with 31 words, some of which are a title, the cell para- C meters and NSYM, the # of symmetry records that follow C 2. NSYM records each defining one equivalent position C 3. a record with 13 words, some of which are the # grid points along C each axis and NORN, defining the orientation of the map C The definition of NORN is : C NORN sections across down C 0 Y X Z C 1 X Y Z C 2 Z X Y C The parameters passed to the main program are IS, IA and ID, giving the C number i of axis XYZ(i) of sections, across and down, respectively; and C NXYZ(i), IORI(i) and IEND(i) , or # grid points and the first and last C point (in grid units) along each axis. C*** C DIMENSION ITLE(20),PARA(6),TRANS(3,24),IROT(2,3,24),MNX(6),NUVW(3) COMMON /GRID/ NXYZ(3), IORI(3), IEND(3), IOUT1(2), IOUT2(2) C READ (MI) ITLE, NOSET, PARA, NSYM, NCENT, LATT, NPIC WRITE (IO,10) (ITLE(I),I=1,16), PARA WRITE (LP,11) ITLE, PARA DO 1000 I=1,NSYM READ (MI) (TRANS(J,I),(IROT(K,J,I),K=1,2),J=1,3) 1000 CONTINUE READ (MI) NUVW,SCALE,NORN,MNX,RSPMX,NBIT WRITE (LP,12) SCALE,MNX,RSPMX OPEN (UNIT=99,FILE='CELL',FORM='FORMATTED',STATUS='NEW') WRITE (99,40) PARA CLOSE (UNIT=99) C C-- Determine the orientation (defined by IS, IA and ID) from NORN : IF (NORN.EQ.0) THEN IS = 2 IA = 1 ID = 3 ELSE IF (NORN.EQ.1) THEN IS = 1 IA = 2 ID = 3 ELSE IF (NORN.EQ.2) THEN IS = 3 IA = 1 ID = 2 ELSE WRITE (IO,30) NORN WRITE (LP,31) NORN STOP ENDIF NXYZ(IS) = NUVW(3) NXYZ(IA) = NUVW(1) NXYZ(ID) = NUVW(2) C C-- Origin and last point on each axis (in grid units) : DO 2000 I=1,3 IORI(I) = 0 IEND(I) = NXYZ(I) - 1 2000 CONTINUE RETURN C 10 FORMAT (T2,'Title read : ',16A4 . /T2,'Cell parameters :',5X,6F8.2) 11 FORMAT (/T40,'TITLE READ: ',20A4 . //T40,'CELL PARAMETERS :',10X,6F8.2) 12 FORMAT (/T40,'SCALE FACTOR :',10X,E10.4 . //T40,'MIN. AND MAX. H,K,L',6X,3(2I5,2X) . //T40,'RESOLUTION :',10X,F8.2) 20 FORMAT (T2,'# grid points along X,Y,Z : ',3I4) 21 FORMAT (/T40,'SCALE FACTOR',10X,F8.4 . //T40,'MIN. AND MAX. INDICES ',10X,6I6 . //T40,'RESOLUTION ',10X,F8.2) 30 FORMAT (//T30,'*** ORIENTATION OF THE OLD MAP CANNOT BE ', . 'DECODED : NORN =',I4,' ***') 31 FORMAT (////T30,'*** ORIENTATION OF THE OLD MAP CANNOT BE ', . 'DECODED : NORN =',I4,' ***') 40 FORMAT ('CELL ',6F10.4) END C C=========================================================================== C SUBROUTINE ORIENT(II,IO,IS,IA,ID,JS,JA,JD,NORIEN,IOU,IOV) C DIMENSION IACR(3), IORIEN(3) COMMON /GRID/ NXYZ(3), IORI(3), IEND(3), IOUT1(2), IOUT2(2) C C-- FORDAP definitions : IORIEN(1) = 100 IORIEN(2) = 10 IORIEN(3) = 1 IACR(1) = 3 IACR(2) = 1 IACR(3) = 1 C C-- New orientation....(no changes to input map allowed!) 1 JS = IS JA = IA JD = ID NORIEN = IORIEN(JS) IF (IACR(JS).NE.JA) NORIEN = -NORIEN IOU = 0 IOV = 0 RETURN C 90 FORMAT (A) END C C============================================================================= C SUBROUTINE RDMFSFOR(MI,IO,IS,IA,ID,RMAX,RMS) PARAMETER (MAPDIM=40000) C COMMON /MAPS/ MAP(MAPDIM) COMMON /GRID/ NXYZ(3), IORI(3), IEND(3), IOUT1(2), IOUT2(2) K1 = 1 KD = NXYZ(IA) - 1 DO 1100 J=1,NXYZ(ID) K2 = K1 + KD READ (MI) (MAP(K),K=K1,K2) DO 1110 K=K1,K2 IF (RMAX.LT.MAP(K)) RMAX = MAP(K) RMS = RMS + MAP(K)*MAP(K) 1110 CONTINUE K1 = K2 + 1 1100 CONTINUE RETURN C 10 FORMAT ('+',T56,'[ reading section #',I3,' ]') 20 FORMAT ('+',T56,24(' ')) END C C============================================================================= C SUBROUTINE SKIP(MI,IO,IS,IA,ID) PARAMETER (MAPDIM=40000) C COMMON /MAPS/ MAP(MAPDIM) COMMON /GRID/ NXYZ(3), IORI(3), IEND(3), IOUT1(2), IOUT2(2) C DO 1100 J=1,NXYZ(ID) READ (MI) 1100 CONTINUE RETURN END C C============================================================================= C SUBROUTINE WRMCTR(MO,IO,LP,OUTFILE,IS,IA,ID,JS,JA,JD,NF,NL, . NORIEN,RMAX,MAPTYP,ISEC) PARAMETER (MAPDIM=40000) C DIMENSION XYZ(3), ISTP(3), MAP(MAPDIM), IOSH(3) COMMON /MAPS/ AMAP(MAPDIM) COMMON /GRID/ NXYZ(3), IORI(3), IEND(3), IOUT1(2), IOUT2(2) COMMON /ORTH/ YESNO EQUIVALENCE (MAP(1),AMAP(1)) CHARACTER*1 XYZ, OUTFILE*14 LOGICAL YESNO DATA XYZ/'X','Y','Z'/ C C-- Some variables written to the CONTUR map : ZERO = 0.0 SEV1 = 1000.0 DUM = 0.0 I1 = 1 C C-- Now calculate the CONTUR step, map limits, no. points : DU = 1.0 / NXYZ(JA) DV = 1.0 / NXYZ(JD) DW = 1.0 / NXYZ(JS) C UMIN = FLOAT(IOUT1(1))/NXYZ(JA) UMAX = FLOAT(IOUT2(1))/NXYZ(JA) VMIN = FLOAT(IOUT1(2))/NXYZ(JD) VMAX = FLOAT(IOUT2(2))/NXYZ(JD) WMIN = FLOAT(NF)/NXYZ(JS) WMAX = FLOAT(NL)/NXYZ(JS) C NU = IOUT2(1)-IOUT1(1)+1 NV = IOUT2(2)-IOUT1(2)+1 NW = NL-NF+1 NMODU = NXYZ(IA) NMODV = NXYZ(ID) IF (ISEC.EQ.NF) WRITE (LP,30) XYZ(JS),NW,WMIN,WMAX,DW, . XYZ(JA),NU,UMIN,UMAX,DU,XYZ(JD),NV,VMIN,VMAX,DV C C-- And here comes the new map : WRITE (IO,10) ISEC W = ISEC*DW WRITE (MO) SEV1, SEV1 WRITE (MO) UMIN,UMAX,DU,VMIN,VMAX,DV,WMIN,WMAX,DW,NV,NU,I1, . NORIEN,DUM,DUM,W K1 = IOUT1(1) K2 = K1+NU-1 DO 1100 JJ=IOUT1(2),IOUT2(2) ! JJ : LINE V = JJ*DV LINE = MOD(JJ,NMODV) K0 = LINE*NXYZ(IA)+1 WRITE (MO) ZERO, ZERO WRITE (MO) V, (FLOAT(MAP(K0+MOD(K,NMODU))),K=K1,K2) C WRITE (3,40) V, (MAP(K0+MOD(K,NMODU)),K=K1,K2) 1100 CONTINUE RETURN C 66 FORMAT (T2,40I3) 10 FORMAT ('+',T56,'[ writing section #',I3,' ]') 30 FORMAT (/T20,'OUTPUT MAP INFO IN FRACTIONAL COORDINATES'// . T58,'LAYOUT',T73,'#POINTS',T94,'FROM ... TO',T114,'STEP SIZE' . //T40,'SECTIONS :',T61,A,T75,I4,T90,F7.4,3X,F7.4,T114,F9.6/ . T40,'ACROSS :',T61,A,T75,I4,T90,F7.4,3X,F7.4,T114,F9.6/T40, . 'DOWN :',T61,A,T75,I4,T90,F7.4,3X,F7.4,T114,F9.6) 40 format (t2,f5.3,':',30i4/6(t7,30i4/t7,30i4)) END C C----------------------------------------------------------------- C SUBROUTINE POSFSF(IUNIT,NW1,ISEC,NV) C ==================================== C C ----------------------------------------------------------- C POSITIONS A FSFOR MAP FILE TO THE BEGINNING OF SECTION ISEC C [12/88 BdV] C ----------------------------------------------------------- C DIMENSION NWDS(27) REWIND IUNIT C C-- SKIP HEADER RECORDS C READ (IUNIT) NWDS,NSYM DO I=1,NSYM+1 READ (IUNIT) ENDDO C C-- SKIP ISEC-NW1 SECTIONS C DO I=1,ISEC-NW1 DO J=1,NV READ (IUNIT,ERR=9999,END=9999) ENDDO ENDDO RETURN C 9999 WRITE (6,10) ISEC STOP 'ERROR POSITIONING MAP FILE' 10 FORMAT (//T20,'*** ERROR TRYING TO FIND MAP SECTION ',I4) END