C- RESECTION CCP4 MAP 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*4 TLINE(256),symtry(512) EQUIVALENCE (TLINE(1),BLINE(1)) real*4 rline(256),TCELL(6),CELL(6),cell2(6) 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),ITRANS(3),nabc(3) EQUIVALENCE (LINT(1),BLINE(1)) CHARACTER*80 TTITLE,TITLE !,TSYMTRY,SYMTRY EQUIVALENCE (TTITLE,RLINE(57)) EQUIVALENCE (TSYMTRY,RLINE(1)) CHARACTER*50 INFILE,OUTFILE,FILENAME LOGICAL*1 INSIDE DIMENSION RMAP(120,120,300) cc******DIMENSION RMAP(NMED,NSLOW,NFAST) FOR INPUT MAP!!! CC******CHANGE 2 PLACES BELOW IN IF BLOCK IF CHANGED HERE! WRITE(6,1001) 1001 FORMAT(//' PROGRAM TO RESECTION CCP4 MAP FILES', @ /' SWAP MEDIUM OR FAST WITH SLOW AXES, WRITE TO NEW.MAP') 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', . access='direct',recl=256) READ(2)RLINE NFAST=LINT(1) NMED=LINT(2) NSLOW=LINT(3) MAPTYPE=LINT(4) DO 290 I=1,6 290 CELL(I)=TCELL(I) NX1=LINT(8) NY1=LINT(9) NZ1=LINT(10) do 292 i=1,3 292 nabc(i)=lint(7+i) write (6,*)'nabc:',(nabc(i),i=1,3) NSYMCHAR=RLINE(55) TITLE=TTITLE iaxfast=lint(17) iaxmed=lint(18) iaxslow=lint(19) dmin=rline(20) dmax=rline(21) dmean=rline(22) nsg=lint(23) nsymchar=lint(24) ddev=RLINE(55) ntitle=lint(56) WRITE(6,*)'NFAST, MEDIUM, SLOW=',NFAST,NMED,NSLOW IF ((NFAST.GT.300).OR.(NMED.GT.120).OR.(NSLOW.GT.120)) THEN STOP 'DIMENSIONED FOR 300,120,120 !' ENDIF WRITE(6,*)'MAPTYPE=',MAPTYPE write(6,*)' IAXFAST, MEDIUM, SLOW =',IAXFAST,IAXMED,IAXSLOW IF WRITE(6,*)'CELL: ',(CELL(I),I=1,6) WRITE(6,*)'NX1, NY1, NZ1=',NX1, NY1, NZ1 WRITE(6,*)'#CHAR IN SYMOPS',NSYMCHAR WRITE(6,*) TITLE READ(2)RLINE c SYMTRY=TSYMTRY c WRITE(6,*)SYMTRY ii=1 do 295 i=1,nsymchar/4 IF (II.GT.256) THEN READ(2,END=1900) RLINE II=1 ENDIF symtry(i)=TLINE(II) !if symtry char*4 c write(6,*)ii,tline(ii) 295 ii=ii+1 c do 296 i=1,nsymchar/80-1 c ii=20*(i-1) c296 write(6,*) ii,(symtry(ii+j),j=1,20) c II=1+NSYMCHAR/4 C !!!!IF NSYMCHAR NOT MULTIPLE OF 4, SHOULD ROUND UP? DO 325 J=1,NSLOW-1 C WRITE(6,*)'SECTION ',J DO 300 I=1,NMED DO 300 K=1,NFAST IF (II.GT.256) THEN READ(2,END=1900) RLINE II=1 ENDIF RMAP(I,J,K)=RLINE(II) C I,J,K CORRESPOND TO NMED, NSLOW, NFAST OR X,Y,Z FOR Y-SECTION MAPS. 300 II=II+1 325 CONTINUE 380 CLOSE(UNIT=2) C NOW WE HAVE THE WHOLE MAP IN MEMORY! C PREPARE SOME THINGS: C MAKE MATRIX FOR NCS SYMMETRY OPERATION ON GRID UNITS. CALL SYMMAT(ARRAY,cell,nabc) OPEN (UNIT=3,FILE='new.map',STATUS='unknown',FORM='UNFORMATTED', & RECL=256, access='direct', CARRIAGECONTROL='NONE') DO 400 I=57,256 400 TLINE(I)=' ' DO 410 I=1,56 410 RLINE(I)=0 LINT(4)=MAPTYPE DO 490 I=1,6 490 TCELL(I)=CELL(I) RLINE(55)=RMYST TTITLE=TITLE rline(20)=DMIN rline(21)=DMAX rline(22)=DMEAN lint(23)=NSG lint(24)=NSYMCHAR RLINE(55)=DDEV lint(56)=NTITLE LINT(8)=NX1 LINT(9)=NY1 LINT(10)=NZ1 IF (ISWAP.EQ.0) THEN C SWAP MED AND SLOW AXES: 500 LINT(1)=NFAST LINT(3)=NMED LINT(2)=NSLOW lint(17)=IAXFAST lint(18)=IAXSLOW lint(19)=IAXMED WRITE(3) RLINE ii=1 DO 600 I=1,nsymchar/4 IF (II.GT.256) THEN write(6,*)' writing a record:' WRITE(3) RLINE II=1 ENDIF TLINE(II)=symtry(I) c write (6,*)ii,tline(ii) 600 ii=ii+1 c II=1+NSYMCHAR/4 DO 700 I=1,NMED C WRITE(6,*) WRITE(6,*)'SECTION ',I DO 700 J=1,NSLOW DO 700 K=1,NFAST IF (II.GT.256) THEN WRITE(3) RLINE II=1 ENDIF c-------------------------------------------------------------- c put into this i,j,kth pixel of new map the product of c density at i,j,kth pixel and density at its symm related c position ii,jj,kk. c call symncs(i,j,k,array,ii,jj,kk) ri(1)=i ri(2)=j ri(3)=k do 652 in=1,3 ro(in)=0 do 650 jn=1,3 650 ro(in)= ro(in)+array(in,jn)*ri(jn) 652 ro(in)= ro(in)+array(in,4) ii=int(ro(1)+0.5) jj=int(ro(2)+0.5) kk=int(ro(3)+0.5) x=(RMAP(I,J,K))*(RMAP(II,JJ,KK)) !!assume average of map is zero, this gives correlation. RLINE(II)=X 700 II=II+1 C------------------------------------------------------------------ WRITE(3) RLINE CLOSE(UNIT=3) STOP 'NORMAL END' 1900 STOP 'END OF FILE BEFORE END OF MAP' END SUBROUTINE SYMMAT(ARRAY,a,b,c,alpha,beta,gamma,na,nb,nc) CHARACTER*40 INFILE,FILENAME LOGICAL EFLG real*8 M1(3,4),M2(3,4),ncs(3,4),cell(3,4),rcell(3,4),pi c REAL*8 P1(3),P2(3),T(3) character*40 string INTEGER*4 MM,NN,I,J pi = 4*atan(1.) c write(6,*)'Enter cell constants a,b,c, alpha, beta, gamma:' c read(5,*)a,b,c, alpha, beta, gamma alpha=alpha*pi/180. beta = beta*pi/180. gamma=gamma*pi/180. cell(1,1)=a cell(1,2)=b*cos(gamma) cell(2,2)=b*sin(gamma) cell(1,3)=c*cos(beta) cell(2,3)=c*(cos(alpha)-cos(beta)*cos(gamma))/sin(gamma) cell(3,3)=sqrt(c*c-(cell(3,1)*cell(3,1)+cell(3,2)*cell(3,2))) do 20 i=1,3 do 20 j=1,3 20 rcell(i,j)=cell(i,j) mm=3 call matinv(rcell,mm,j) c read operator into m1 5 WRITE(6,*) 'ENTER FILENAME rotation-trans op (.o) on orthog Angst:' READ(5,1200) INFILE 1200 FORMAT(A) IF (INFILE(:5).EQ.' ') STOP INQUIRE(FILE=INFILE, NAME=FILENAME,EXIST=EFLG) IF (.NOT.EFLG) GOTO 5 OPEN (UNIT=2,FILE=INFILE,STATUS='OLD',READONLY) c READ(2,*) '.SPACE_GROUP_operators r 12 (3f12.7)' READ(2,'(A40)') string do 1150 j=1,4 1150 read(2,*)(ncs(I,J),I=1,3) 1151 FORMAT(3f14.7) CLOSE(UNIT=2) MM=3 NN=4 CALL MATPRN (ncs,MM,NN) c multiply the three matrices together: ncs(i,j) is noncryst sym operator; j=4 => translation. do 30 i=1,3 do 30 j=1,3 m1(i,j)=0 do 30 k=1,3 30 m1(i,j)=m1(i,j)+ncs(i,k)*cell(k,j) do 40 i=1,3 do 40 j=1,3 m2(i,j)=0 do 40 k=1,3 40 m2(i,j)=m2(i,j)+rcell(i,k)*m1(k,j) do 50 i=1,3 m2(i,4)=0 do 50 k=1,3 50 m2(i,4)=m2(i,4)+rcell(i,k)*ncs(k,4) WRITE(6,*)' NCS OPERATOR (fractional):' WRITE(6,*)' ROTATE MATRIX' DO 1100 I=1,3 1100 WRITE(6,9) (M2(I,J),J=1,3) WRITE(6,*)' TRANSLATE ' WRITE(6,9) (M2(I,4),I=1,3) 9 FORMAT(' ',3F14.6) c multiply by grid on right and 1/grid on left to work directly in grid units: do 1210 i=1,3 do 1200 j=1,3 1200 m2(i,j)=m2(i,j)*nabc(j)/nabc(i) 1210 m2(i,4)=m2(i,4)/nabc(i) WRITE(6,*)' NCS OPERATOR (grid units):' WRITE(6,*)' ROTATE MATRIX' DO 1250 I=1,3 1250 WRITE(6,9) (M2(I,J),J=1,3) WRITE(6,*)' TRANSLATE ' WRITE(6,9) (M2(I,4),I=1,3) return end SUBROUTINE MATPRN(ARRAY,M,N) REAL*8 ARRAY(3,4) WRITE(6,*)'M, N=',M,N DO 10 I=1,M 10 WRITE (6,50) (ARRAY(I,J),J=1,N) 50 format(8f10.5) RETURN END SUBROUTINE MATINV (ARRAY, NORDER, DET) DOUBLE PRECISION ARRAY, AMAX, SAVE 2 DIMENSION ARRAY(3,4), IK (10), JK (10) 10 DET = 1 11 DO 100 K=1, NORDER C C FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX C AMAX=0 21 DO 30 I=K, NORDER DO 30 J=K, NORDER 23 IF (DABS (AMAX) - DABS (ARRAY(I,J) ) ) 24, 24, 30 24 AMAX=ARRAY(I,J) IK(K) =I JK(K) =J 30 CONTINUE C C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K) C 31 IF (AMAX) 41, 32, 41 32 DET=0. GO TO 140 41 I=IK (K) IF (I-K) 21, 51, 43 43 DO 50 J=1, NORDER SAVE= ARRAY(K,J) ARRAY(K,J)= ARRAY(I,J) 50 ARRAY(I,J)= -SAVE 51 J=JK (K) IF (J-K) 21, 61, 53 53 DO 60 I=1, NORDER SAVE= ARRAY(I,K) ARRAY(I,K)= ARRAY(I,J) 60 ARRAY(I,J)= -SAVE C C ACCUMULATE ELEMENTS OF INVERSE MATRIX C 61 DO 70 I=1, NORDER IF(I-K) 63, 70, 63 63 ARRAY(I,K) = -ARRAY(I,K) / AMAX 70 CONTINUE 71 DO 80 I=1, NORDER DO 80 J=1, NORDER IF (I-K) 74, 80, 74 74 IF (J-K) 75, 80, 75 75 ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K) *ARRAY(K,J) 80 CONTINUE 81 DO 90 J=1, NORDER IF (J-K) 83, 90, 83 83 ARRAY(K,J) = ARRAY(K,J) / AMAX 90 CONTINUE ARRAY(K,K) = 1. / AMAX 100 DET=DET *AMAX C C RESTORE ORDERING OF MATRIX C 101 DO 130 L=1, NORDER K= NORDER - L +1 J= IK (K) IF (J-K) 111, 111, 105 105 DO 110 I=1, NORDER SAVE= ARRAY(I,K) ARRAY(I,K)= -ARRAY(I,J) 110 ARRAY(I,J)= SAVE 111 I= JK(K) IF (I-K) 130, 130, 113 113 DO 120 J=1, NORDER SAVE= ARRAY(K,J) ARRAY(K,J) = -ARRAY(I,J) 120 ARRAY(I,J)= SAVE 130 CONTINUE 140 RETURN END