C index in rline C V -this is taken care of by do loop structure of ccpmap3 C coord in output map C V -subt 1 from each index so first is at 0,0,0 C divide by extent and multiply by cell edge for AU C angst pos in output map C V -mult by matrix and translate C angst pos in input map C V - divide by cell edge and multiply by extent C coord in input map C V - interplate between pixels for density? C Take density at this location and assign it to the curr pos in new map C C From old map file get nfast, med, slow; cell param C need same for new map from input stream C- READ 1024 BYTES AT A TIME. C- READ INTO ARRAY COMMON WITH 256 REAL*4, 256 INT*4, 1 BYTE STRINGS C- PICK OUT THE VALUES FOR A CCP4 DENSITY MAP, BE THEY CHARACTER, INTEGER, C- OR REAL (DENSITIES ARE REAL). CHARACTER*1 BLINE(1024) CHARACTER*36 RECLINE(30) EQUIVALENCE (RECLINE(1),BLINE(1)) real*4 rline(256),TCELL(6),CELL1(6),cell2(6),TRANS(3) REAL*4 RMX(3,3),RMY(3,3),RMZ(3,3),TEMP(3,3),PROD(3,3),X1(3),X2(3) EQUIVALENCE (RLINE(1),BLINE(1)),(TCELL(1),RLINE(11)) INTEGER*4 LINT(256) EQUIVALENCE (LINT(1),BLINE(1)) CHARACTER*80 TTITLE,TITLE,TITLE2,TSYMTRY,SYMTRY EQUIVALENCE (TTITLE,RLINE(56)) EQUIVALENCE (TSYMTRY,RLINE(1)) CHARACTER*50 INFILE,OUTFILE,FILENAME LOGICAL*1 INSIDE DIMENSION RMAP(100,100,300) WRITE(6,1001) 1001 FORMAT(//' PROGRAM TO READ CCP4 MAP FILES') 5 WRITE(6,1002) 1002 FORMAT(/'$Input filename (CR to quit): ') READ(5,1003) INFILE 1003 FORMAT(A) C IF (INFILE(:5).EQ.' ') STOP C INQUIRE(FILE=INFILE, NAME=FILENAME,EXIST=EFLG) C IF (.NOT.EFLG) GOTO 5 OPEN (UNIT=2,FILE=INFILE,STATUS='OLD',FORM='UNFORMATTED') READ(2)(RLINE(I),I=1,256) NFAST1=LINT(1) NMED1=LINT(2) NSLOW1=LINT(3) maptype=LINT(4) nsfast=lint(5) nsmed=lint(6) nsslow=lint(7) DO 290 I=1,6 290 CELL1(I)=TCELL(I) C These dimensions are sampling dimensions, corresp to full unit cell: NX1=LINT(8) NY1=LINT(9) NZ1=LINT(10) iaxfast=lint(17) iaxmed=lint(18) iaxslow=lint(19) dmin=rline(20) dmax=rline(21) dmean=rline(22) nsg1=lint(23) nmyst=lint(24) ddev=RLINE(55) ntitle=lint(56) TITLE=TTITLE WRITE(6,*)'NFAST1, MEDIUM1, SLOW1=',NFAST1,NMED1,NSLOW1 WRITE(6,*)'MAP TYPE=',MAPTYPE WRITE(6,*)'START (FAST, MED, SLOW)=',NSFAST,NSMED,NSSLOW WRITE(6,*)'CELL: ',(CELL1(I),I=1,6) WRITE(6,*)'SAMPLING(NX, NY, NZ)=',NX1, NY1, NZ1 WRITE(6,*)'AXIS PERMUTAION (FAST, SLOW, MEDIUM)=',IAXFAST, & IAXMED,IAXSLOW WRITE(6,*)'DMIN, DMAX, DMEAN, DDEV=',DMIN,DMAX,DMEAN,DDEV WRITE(6,*)'SPACE GROUP #',NSG1 WRITE(6,*)'NUMBER OF TITLES:',NTITLE WRITE(6,*) TITLE READ(2)(RLINE(I),I=1,256) SYMTRY=TSYMTRY WRITE(6,*)SYMTRY II=21 DO 325 J=1,NSLOW1 C WRITE(6,*)'SECTION ',J DO 300 I=1,NMED1 DO 300 K=1,NFAST1 IF (II.GT.256) THEN READ(2,END=1900)(RLINE(IJ),IJ=1,256) II=1 ENDIF RMAP(I,J,K)=RLINE(II) C I,J,K CORRESPOND TO X,Y,Z, OR NMED, NSLOW, NFAST 300 II=II+1 325 CONTINUE 380 CLOSE(UNIT=2) C NOW WE HAVE THE WHOLE MAP IN MEMORY! C PREPARE SOME THINGS: write(6,*)'space group number of new cell?' read(5,*) nsg2 WRITE(6,*)'ENTER CELL PARAMETERS FOR NEW CELL:' READ(5,*) (CELL2(I),I=1,6) WRITE(6,*)'ENTER NUMBER OF PIXELS ALONG A,B,C AXES' READ(5,*)NX2,NY2,NZ2 NFAST2=NZ2 NMED2=NX2 NSLOW2=NY2 NORIEN=2 pi=4*atan(1.) C 'read rot angles: WRITE(6,*)'OPERATION TO TAKE NEW CELL TO OLD, i.e. INVERSE:' WRITE(6,*)'ENTER ANGLES FOR ROTATION ABOUT X,Y,Z (STANDARD ORDER):' READ(5,*) rotx,roty,rotz rotx=rotx*pi/180 roty=roty*pi/180 rotz=rotz*pi/180 C READ TRANSLATION WRITE(6,*)'ENTER X,Y,Z ELEMENTS OF TRANSLATION ', @ '(AU) TO BE APPLIED IN NEW CELL AFTER ROTATION:' READ(5,*) (TRANS(I),I=1,3) C SET UP ROTATION MATRIX: rmx(1,1)=1 rmx(2,2)=COS(rotx) rmx(3,2)=SIN(rotx) rmx(2,3)=-rmx(3,2) rmx(3,3)=rmx(2,2) rmy(2,2)=1 rmy(1,1)=COS(roty) rmy(3,1)=SIN(roty) rmy(1,3)=-rmy(3,1) rmy(3,3)=rmy(1,1) rmz(3,3)=1 rmz(1,1)=COS(rotz) rmz(2,1)=SIN(rotz) rmz(1,2)=-rmz(2,1) rmz(2,2)=rmz(1,1) c multiply the matrices together, put result in prod(). do 390 i=1, 3 do 390 K=1, 3 temp(i,K)=0 do 390 j=1, 3 390 temp(i,K)=temp(i,K)+rmy(i,j)*rmz(j,K) do 394 i=1, 3 do 394 K=1, 3 prod(i,K)=0 do 394 j=1, 3 394 prod(i,K)=prod(i,K)+rmx(i,j)*temp(j,K) OPEN (UNIT=3,FILE='NEW.MAP',STATUS='NEW',FORM='UNFORMATTED', & RECL=256, RECORDTYPE='FIXED', CARRIAGECONTROL='NONE') DO 400 I=1,256 400 RLINE(I)=0 LINT(1)=NFAST2 LINT(2)=NMED2 LINT(3)=NSLOW2 LINT(4)=maptype lint(5)=nsfast lint(6)=nsmed lint(7)=nsslow LINT(8)=NX2 LINT(9)=NY2 LINT(10)=NZ2 DO 490 I=1,6 490 TCELL(I)=CELL2(I) lint(17)=iaxfast lint(18)=iaxmed lint(19)=iaxslow rline(20)=dmin rline(21)=dmax rline(22)=dmean lint(23)=nsg2 lint(24)=nmyst rline(55)=ddev lint(56)=ntitle TTITLE=TITLE2 WRITE(3)(RLINE(I),I=1,256) WRITE(6,*)'NFAST, MEDIUM, SLOW=',NFAST2,NMED2,NSLOW2 WRITE(6,*)'maptype=',maptype WRITE(6,*)'CELL: ',(CELL2(I),I=1,6) WRITE(6,*)'NX2, NY2, NZ2=',NX2, NY2, NZ2 WRITE(6,*)'note dmin, dmax, dmean, and ddev are not valid!' WRITE(6,*) TITLE2 DO 493 I=1,20 493 RLINE(I)=0 TSYMTRY=SYMTRY II=21 DO 1500 J=1,NSLOW2 C WRITE(6,*) WRITE(6,*)'SECTION ',J DO 1500 I=1,NMED2 DO 1500 K=1,NFAST2 IF (II.GT.256) THEN WRITE(3)(RLINE(IJ),IJ=1,256) II=1 ENDIF C GET A VALUE FOR THE NEXT (I,J,K'TH) PIXEL C WRITE(6,*)'I,J,K-',I,J,K C CALCULATE X2()=POSITION IN ANGST IN NEW CELL X2(1)=CELL2(1)*(I-1)/NX2 X2(2)=CELL2(2)*(J-1)/NY2 X2(3)=CELL2(3)*(K-1)/NZ2 C WRITE(6,*)'X2=',(X2(I5),I5=1,3),' AU ' C OPERATE ON X2 WITH TRANSLATION TO GET POSITION BEFORE TRANSLATION DO 1000 IJ=1,3 1000 X2(IJ)=X2(IJ)-TRANS(IJ) C WRITE(6,*)'X2=',(X2(I5),I5=1,3) C OPERATE ON X2 WITH MATRIX TO GET POSITION BEFORE ROTATING -> X1() DO 1100 IJ=1,3 X1(IJ)=0 DO 1100 JJ=1,3 1100 X1(IJ)=X1(IJ)+PROD(IJ,JJ)*X2(JJ) C WRITE(6,*)'X1=',(X1(I5),I5=1,3),' AU IN OLD CELL' C TRANSLATE BY CELL EDGE TO BRING INTO CELL IF OUTSIDE: C DO 1200 IJ=1,3 C IF (X1(IJ).LT.0) X1(IJ)=X1(IJ)+CELL1(IJ) C1200 IF (X1(IJ).GT.CELL1(IJ)) X1(IJ)=X1(IJ)-CELL1(IJ) C WRITE(6,*)'X1=',(X1(I5),I5=1,3),' INSIDE CELL' C CHECK IF INSIDE CELL, SET ZERO IF OUT: INSIDE=.TRUE. DO 1200 IJ=1,3 IF (X1(IJ).LT.0) INSIDE=.FALSE. 1200 IF (X1(IJ).GT.CELL1(IJ)) INSIDE=.FALSE. IF (.NOT. INSIDE) THEN RLINE(II)=0 GOTO 1500 ENDIF C CONVERT TO PIXEL UNITS X1(1)=X1(1)*NX1/CELL1(1) X1(2)=X1(2)*NY1/CELL1(2) X1(3)=X1(3)*NZ1/CELL1(3) C WRITE(6,*)'X1=',(X1(I5),I5=1,3),' PIXELS' C ROUND OFF TO WHOLE PIXEL (MAYBE INTERPOLATE BETWEEN PIXELS?) I1=1+INT(X1(1)+.5) J1=1+INT(X1(2)+.5) K1=1+INT(X1(3)+.5) C WRITE(6,*)'I1,J1,K1=',I1,J1,K1 RLINE(II)=RMAP(I1,J1,K1) 1500 II=II+1 WRITE(3)(RLINE(I),I=1,256) STOP 'NORMAL END' 1900 STOP 'END OF FILE BEFORE END OF MAP' END C STILL NEED TO APPLY SYMETRY OPERATIONS! USE ANOTHER PROGRAM, READ IN THE MAP ALL INTO MEMORY, C STEP THROUGH THE PIXELS OF ONE ASU. FOR EACH PIXEL, EXAMINE ALL EQUIVLENT PIXELS. SET ALL TO VALUE OF DENSEST. C ONLY ONE OF THE EQUIV PIXELS SHOULD HAVE ANY DENSITY BEFORE THIS?