c read PDB file in orthogonal system (later read any pdb, convert to orth axes) c no, read RDI. convert pdb to rdi w pdb2rdi AFTER REMOVING nonatom records c operate with specified rotation + translation c put into arbitrary unit cell (given orth axis relative to cell vectors) c (write out this file for use by sfclc) c read MEP file c check for overlap of protein between asu's: c for each atom c operate on atom with each MEP operation c check atom in new position for overlap with each atom of primary position c (If primary position does not overlap with any other, that ensures c there will be no overlap between an positions?) REAL*4 V(3),W(3),RM(3,3),A,B,ANGLE,COSA,SINA,SYMOP(3,4,12) REAL*4 CELLA,CELLB,CELLC,ALPHA,BETA,GAMMA,CELL(3,3) INTEGER I,J,K,I1,I2,I3 CHARACTER*15 RESI CHARACTER*40 FNAME CHARACTER*1 AXIS character*80 astring WRITE(6,*) 'PROGRAM TO CHECK .RDI FILE FOR COLLISION BETWEEN' WRITE(6,*) 'ASYMMETRIC UNITS. .RDI FILE HAS ORTH COORDINATES' WRITE(6,*) 'WRITE NEW .RDI FILE WITH COOR AU ALONG CELL AXES' 30 TYPE *,'ENTER FILE NAME OF RDI FILE TO CHECK:' 35 read(5,40)FNAME IF (FNAME.EQ.'Q') GOTO 501 40 FORMAT (A40) OPEN(UNIT=2,FILE=FNAME,STATUS='OLD') ' TYPE *,'ENTER FILE NAME FOR OUTPUT:' ' read(5,40)FNAME ' OPEN(UNIT=3,FILE=FNAME,STATUS='NEW') WRITE(6,*)'ENTER A B C ALPHA BETA GAMMA:' READ(5,*) CELLA,CELLB,CELLC,ALPHA,BETA,GAMMA CELL(1,1)=CELLA CELL(2,1)=0 CELL(3,1)=0 CELL(1,2)=CELLB*COS(GAMMA) CELL(2,2)=CELLB*SIN(GAMMA) CELL(3,2)=0 CELL(1,3)=CELLC*COS(BETA) CELL(2,3)=CELLC*COS(ALPHA)*SIN(BETA) CELL(3,3)=CELLC*SIN(ALPHA)*SIN(BETA) DO 405 I=1,3 DO 405 J=1,3 405 RCELL(I,J)=CELL(I,J) N=3 CALL MATINV(RCELL,N,DET) C IF V IS COOORD IN FRACTION ALONG UNIT CELL => X = CELL * V C SO V CAN BE OBTAINED AS V = RCELL * X C MEP operators work on V, not on X. WRITE(6,*) ENTER NAME OF .MEP FILE WITH SYM OPS READ(5,40)FNAME OPEN (UNIT=4,FILE=FNAME,STATUS='OLD') READ(4,*)NSYM DO 41 K=1,NSYM 41 READ(4,*) ((SYMOP(I,J,K),I=1,3),J=1,4) CLOSE (UNIT=4) WRITE(6,*)'ENTER AXIS TO ROTATE ABOUT (X, Y, Z) AND ANGLE . FOR ROTATION' READ(5,42)AXIS,ANGLE 42 FORMAT(A1,F15.3) PI=4*ATAN(1.) WRITE(6,*)'PI=',PI,' ANGLE=',ANGLE COSA=COS(ANGLE*PI/180.) SINA=SIN(ANGLE*PI/180.) DO 45 I=1,3 DO 45 J=1,3 RX(I,J)=0 RY(I,J)=0 45 RZ(I,J)=0 RX(1,1)=1 RX(2,2)=COSA RX(3,3)=COSA RX(3,2)=SINA RX(2,3)=-SINA RY(2,2)=1 RY(1,1)=COSA RY(3,3)=COSA RY(3,1)=SINA RY(1,3)=-SINA RZ(3,3)=1 RZ(1,1)=COSA RZ(2,2)=COSA RZ(2,1)=SINA RZ(1,2)=-SINA DO 48 I=1,3 48 WRITE(6,*)(RM(I,J),J=1,3) 49 format (a80) read(2,49) astring write(3,*) astring read(2,49) astring write(6,*)'title from input file:',astring write(3,*) 'rotated about axis ',axis,' by ',angle read(2,49) astring write(3,*) astring C 20.000 32.000 4.000 0.000 2 1 0 1.00 ALA 1 N 50 format (4f10.3,3I5,f8.2,A15) C50 format (4f10.3,3I5,f8.2,' ALA 1 N ') C50 format ('HETATM',I5,' N ALA A 1 ',3f8.3,2f6.2) 55 read(2,50,END=501) (v(i),i=1,3),a,I1,I2,I3,B,RESI write (6,*) (v(i),i=1,3) c operate on atom, now in angstroms along orthogonal coord, w rot,trans do 90 i=1,3 w(i)=rm(i,4) do 90 j=1,3 90 w(i)=w(i)+rm(i,j)*v(j) write (6,*) (w(i),i=1,3) c now operate on repositioned molecule to express its position in angstroms c along unit cell (Faster to form prod of 2 matrices and mult each atom c only once, but keeping separate should help debugging.) C IF V IS COOORD IN FRACTION ALONG UNIT CELL => X = CELL * V C SO V CAN BE OBTAINED AS V = RCELL * X C MEP operators work on V, not on X. do 92 i=1,3 v(i)=0 do 92 j=1,3 92 v(i)=v(i)+rcell(i,j)*w(j) write (6,*) (v(i),i=1,3) c write it back out for further use: cC write(3,50) i1,i2,(w(i),i=1,3),a,b c WRITE(3,50) (W(i),i=1,3),a,I1,I2,I3,B,RESI c GOTO 55 501 CONTINUE END c operate on atom with symmetry operator: for ksym=2,nsym do 90 i=1,3 w(i)=symop(i,4,ksym) do 90 j=1,3 90 w(i)=w(i)+rm(i,j,ksym)*v(j) SUBROUTINE MATPRN(ARRAY,M,N) C DOUBLE PRECISION ARRAY REAL*8 ARRAY(10,10) DO 10 I=1,M 10 WRITE (6,*) (ARRAY(I,J),J=1,N) RETURN END SUBROUTINE MATINV (ARRAY, NORDER, DET) DOUBLE PRECISION ARRAY, AMAX, SAVE 2 DIMENSION ARRAY(10,10), 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