PROGRAM PREPHSP PARAMETER (MAPDIM=50000) C** C This program reads a difference Patterson, and writes out a map in C FFT-format. This map can then be input to HASSP. Currently accepted C map types are stored in array MAPNAMES below. C C HASSP requires the map to be oriented as z sections, with x across C and y down the page. If the input map is oriented differently, this C program will stop. The accepted orientation codes for each map type C are stored in array NORC - compare the comments at the start of the C subroutines that read the header records (SUBROUTINE HDRxxxx). C C The name and size of the map file to be prepared for HASSP is read C from file HASSP.IN. For completeness sake a description of the first C three records, which are read by this program ( subroutine GETINP ), C follows (see also the HASSP manual) : C C Rec. #1 : NAME OF THE FILE CONTAINING MATRICES OF EQUIV. POSITIONS C #2 : NAME OF THE MAP FILE TO BE READ BY HASSP C #3 : GRID(1), GRID(2), GRID(3), NPATXE, NPATYE, NPATZE C ( # grid points along x,y,z; ending grid points for map) C** C DIMENSION MAPNAMES(5), NORC(5), MAP(MAPDIM) COMMON /MAPS/ AMAP(MAPDIM) C-- Variables describing the input map : COMMON /GRID/ NGRID(3), NXYZ(3), IORI(3), IEND(3) EQUIVALENCE (MAP(1),AMAP(1)) CHARACTER*40 FILOUT, FILIN CHARACTER*10 MAPNAMES DATA II/5/, IO/6/, LP/3/, MI/1/, MO/2/ C C-- NOPT is the number of recognized options. MAPNAMES(NOPT) con- C tains the names of the allowed input map types. NORC(NOPT) C contains the acceptable orientation code for each map type : DATA NOPT/2/, MAPNAMES /' FSFOR.MAP','CONTUR.MAP',' ',' ',' '/ DATA NORC/2,1,0,0,0/ C C---------------------------------------------------------------- C C-- Open input map file (the user can choose between the maps C indicated in array MAPNAMES) : C WRITE (IO,5) 2 DO 1000 I=1,NOPT WRITE (IO,10) MAPNAMES(I), I 1000 CONTINUE WRITE (IO,15) READ (II,*) MAPTYP IF (MAPTYP.LT.1 .OR. MAPTYP.GT.NOPT) GOTO 2 1 WRITE (IO,20) MAPNAMES(MAPTYP) READ (II,30) FILIN IF (FILIN.EQ.' ') FILIN = MAPNAMES(MAPTYP) OPEN (UNIT=MI,FILE=FILIN,STATUS='OLD',FORM='UNFORMATTED', . READONLY,ERR=999) C C-- Read the input file to determine the size of the output map, C-- then read the header records from the input map : C CALL GETINP(IO,NXH,NYH,NZH,NXNW,NYNW,NZNW,FILOUT) GOTO (100,110) MAPTYP C-- FSFOR map : 100 CALL HDRFSFOR(IO,LP,MI,IS,IA,ID,NORN) GOTO 800 C-- CONTUR map : 110 CALL HDRCTR(IO,LP,MI,IS,IA,ID,NORN) REWIND MI GOTO 800 C 800 WRITE (IO,35) NXYZ C C-- 1. Is orientation x across, y down, z sections ? C-- 2. Is the old grid fine enough ? C-- 3. Is the origin included ? C-- 4. Is the old grid consistent with the specification in HASSP.IN ? C IF (NORN.NE.NORC(MAPTYP)) THEN WRITE (IO,40) NORN STOP 'Orientation not z sections, x across and y down' ENDIF IF (NXYZ(1).LT.NXNW.OR.NXYZ(2).LT.NYNW.OR.NXYZ(3).LT.NZNW) THEN WRITE (IO,50) NXYZ, NXNW, NYNW, NZNW STOP 'Input map too small' ENDIF IF (IORI(1).NE.0 .OR. IORI(2).NE.0 .OR. IORI(3).NE.0) THEN WRITE (IO,55) IORI STOP 'Origin not included' ENDIF IF (NGRID(1).NE.NXH .OR. NGRID(2).NE.NYH .OR. NGRID(3).NE.NZH) . WRITE (IO,60) NGRID, NXH, NYH, NZH C C-- Open the output map file. Read the required sections C-- from the input map, and write them to the new map file : C WRITE (IO,80) FILOUT, NXNW, NYNW, NZNW OPEN (UNIT=MO,FILE=FILOUT,STATUS='NEW',FORM='UNFORMATTED') C DO 2000 I=1,NZNW GOTO (200,210) MAPTYP C-- FSFOR map : 200 CALL RDMFSFOR(MI,IO,IS,IA,ID) GOTO 900 C-- CONTUR map : 210 CALL RDMCTR(MI,IO,LP,IS,IA,ID) GOTO 900 C-- Write out this section : 900 CALL WRNEWMAP(MO,NXNW,NYNW,IS,IA,ID,NXYZ(IA),MAPTYP) 2000 CONTINUE C C C-- That 's all ! C STOP 'Normal end' C------------------------------ C C-- Error while trying to open the input file : C 999 WRITE (IO,70) FILIN GOTO 1 C 5 FORMAT (/T2,'What type of input map do you have ?') 10 FORMAT (T20,A,' : MAPTYP =',I2) 15 FORMAT (T2,'Enter MAPTYP',T52,': ',$) 20 FORMAT (/T2,'Enter the name of the input map file', . T40,'[',A,'] : ',$) 30 FORMAT (A) 35 FORMAT (T2,'#points along x,y,z on your old map file:',3I6) 40 FORMAT (/T10,'*** HASSP needs z-sections ; NORN =',I5,' ***'/) 50 FORMAT (/T10,'*** Old map has fewer grid points than required ***' . /T10,'*** old : ',3I6,' required :',3I6,'***'/) 55 FORMAT (/T10,'*** Origin is not at 0 0 0 but at ',3I3,' ***'/) 60 FORMAT (/T2,'!!! WARNING : the grid specified in HASSP.IN is', . ' wrong !!!'/T16,'NX, NY, NZ from the FSFOR map :',3I6 . /T16,'NX, NY, NZ from HASSP.IN :',3I6/) 70 FORMAT (T2,'*** Error while trying to open ',A) 80 FORMAT (T2,'New map file is ',A/T2,'with',I4,' grid points along', . ' x,',I4,' along y and',I4,' along z') END C C=========================================================================== C SUBROUTINE HDRFSFOR(IO,LP,MI,IS,IA,ID,NORN) C 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 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 - NGRID : #grid points per unit cell (along crystallographic x,y,z) C - NXYZ : #grid points actually on the map (crystallographic x,y,z) C - IORI : first grid point on the map (crystallographic x,y,z) C - IEND : last grid point on the map (crystallographic x,y,z) C C*** C DIMENSION ITLE(20),PARA(6),TRANS(3,24),IROT(2,3,24),MNX(6),NUVW(3) COMMON /GRID/ NGRID(3), NXYZ(3), IORI(3), IEND(3) C READ (MI) ITLE, NOSET, PARA, NSYM, NCENT, LATT, NPIC WRITE (IO,10) (ITLE(I),I=1,16), 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 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 STOP ENDIF NGRID(IS) = NUVW(3) NGRID(IA) = NUVW(1) NGRID(ID) = NUVW(2) NXYZ(IS) = NGRID(IS) NXYZ(IA) = NGRID(IA) NXYZ(ID) = NGRID(ID) 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,'Cellparameters :',6F8.2) 30 FORMAT (//T30,'*** ORIENTATION OF THE OLD MAP CANNOT BE ', . 'DECODED : NORN =',I4,' ***') END C C=========================================================================== C SUBROUTINE HDRCTR(IO,LP,MI,IS,IA,ID,NORIEN) C*** C Reads the header records from a FORDAP map file : C 1. a record containing (two) dummy values C 2. a record with the map limits and orientation (16 words). C The orientation of the old map is defined by NORIEN : C NORIEN sections across down C 1 Z X Y C 10 Y X Z C 100 X Z Y C C Passed to the main program are IS, IA and ID, giving the number i in C XYZ(i) of the sectioning, across and down axis , respectively ; and C - NGRID : #grid points per unit cell (along crystallographic x,y,z) C - NXYZ : #grid points actually on the map (crystallographic x,y,z) C - IORI : first grid point on the map (crystallographic x,y,z) C - IEND : last grid point on the map (crystallographic x,y,z) C C*** C DIMENSION IORIEN(3), IACR(3) COMMON /GRID/ NGRID(3), NXYZ(3), IORI(3), IEND(3) C C-- FORDAP definitions : IORIEN(1) = 100 IORIEN(2) = 10 IORIEN(3) = 1 IACR(1) = 3 IACR(2) = 1 IACR(3) = 1 C C-- Read dummy record and map limits : READ (MI) READ (MI) UMIN,UMAX,DU,VMIN,VMAX,DV,WMIN,WMAX,DW,NV,NU,NW, . NORIEN,DUM1,DUM2,W C C-- Determine orientation of old map : IS = 0 DO 1000 I=1,3 IF (ABS(NORIEN).EQ.IORIEN(I)) IS = I 1000 CONTINUE IF (IS.LT.1) THEN WRITE (IO,10) NORIEN WRITE (IO,11) NORIEN STOP ENDIF IA = IACR(IS) IF (NORIEN.LT.0) IA = 6 - IS - IA ID = 6 - IS - IA C C-- Calculate limits in grid units : NGRID(IA) = NINT ( 1.0 / DU ) NGRID(ID) = NINT ( 1.0 / DV ) NGRID(IS) = NINT ( 1.0 / DW ) IORI(IA) = NINT ( UMIN * NGRID(IA) ) IORI(ID) = NINT ( VMIN * NGRID(ID) ) IORI(IS) = NINT ( WMIN * NGRID(IS) ) IEND(IA) = NINT ( UMAX * NGRID(IA) ) IEND(ID) = NINT ( VMAX * NGRID(ID) ) IEND(IS) = NINT ( WMAX * NGRID(IS) ) NXYZ(IS) = ABS ( IORI(IS)-IEND(IS) ) + 1 NXYZ(IA) = ABS ( IORI(IA)-IEND(IA) ) + 1 NXYZ(ID) = ABS ( IORI(ID)-IEND(ID) ) + 1 C 10 FORMAT(//T30,'*** ORIENTATION CANNOT BE DETERMINED - NORIEN =',I8) 11 FORMAT(////T30,'*** ORIENTATION CANNOT BE DETERMINED', . ' - NORIEN =',I8,' ***') RETURN END C C=========================================================================== C SUBROUTINE RDMFSFOR(MI,IO,IS,IA,ID) PARAMETER (MAPDIM=50000) C COMMON /MAPS/ MAP(MAPDIM) COMMON /GRID/ NGRID(3), NXYZ(3), IORI(3), IEND(3) K1 = 1 KD = NXYZ(IA) - 1 DO 1100 J=1,NXYZ(ID) K2 = K1 + KD READ (MI) (MAP(K),K=K1,K2) K1 = K2 + 1 1100 CONTINUE RETURN C 10 FORMAT ('+',T56,'[ reading section #',I3,' ]') 20 FORMAT ('+',T56,24(' ')) END C C============================================================================= C SUBROUTINE RDMCTR(MI,IO,LP,IS,IA,ID) PARAMETER (MAPDIM=50000) C COMMON /MAPS/ AMAP(MAPDIM) COMMON /GRID/ NGRID(3), NXYZ(3), IORI(3), IEND(3) C NPNTS = NXYZ(IA) NLIN = NXYZ(ID) NSEC = NXYZ(IS) K1 = 1 KD = NPNTS - 1 C C-- Read a section (it has lots of dummy records !) : READ (MI) READ (MI) DO 1000 J=1,NLIN K2 = K1 + KD READ (MI) READ (MI) V, (AMAP(K),K=K1,K2) K1 = K2 + 1 1000 CONTINUE C RETURN C 10 FORMAT ('+',T56,'[ reading section #',I3,' ]') END C C============================================================================= C SUBROUTINE WRNEWMAP(MO,NXNW,NYNW,IS,IA,ID,NXT,MAPTYP) PARAMETER (MAPDIM=50000) C DIMENSION MAP(MAPDIM) COMMON /MAPS/ AMAP(MAPDIM) EQUIVALENCE (MAP(1),AMAP(1)) C C-- NXT : total # points across (x); only NXNW points per line C are written to the new map; only NYNW lines per section : K1 = 1 KD = NXNW - 1 DO 1100 J=1,NYNW K2 = K1 + KD IF (MAPTYP.EQ.1) THEN WRITE (MO) (FLOAT(MAP(K)),K=K1,K2) ELSE WRITE (MO) (AMAP(K),K=K1,K2) ENDIF K1 = K1 + NXT 1100 CONTINUE RETURN END C C========================================================================== C SUBROUTINE GETINP(IO,NXH,NYH,NZH,NXNW,NYNW,NZNW,FILOUT) CHARACTER*(*) FILOUT C OPEN (UNIT=2,FILE='HASSP.IN',STATUS='OLD',FORM='FORMATTED', . READONLY, ERR=999) READ (2,10) READ (2,20) FILOUT READ (2,10) NXH, NYH, NZH, NXNW, NYNW, NZNW C C-- Add 1 to get the total # points in each direction C (since 0,0,0 is always included in the new map!): NXNW = NXNW + 1 NYNW = NYNW + 1 NZNW = NZNW + 1 C RETURN C 999 WRITE (IO,30) STOP 'HASSP.IN not found' C 10 FORMAT (6I5) 20 FORMAT (A) 30 FORMAT (/T10,'*** Error while trying to open HASSP.IN ***') END