CHARACTER*40 INFILE,FILENAME LOGICAL EFLG real*8 M1(3,4),M2(3,4),M3(3,4),M(3,4),ACC(3,4) REAL*8 P1(3),P2(3),T(3) character*40 string INTEGER*4 MM,NN,I,J c read operator into m1 5 WRITE(6,*) 'ENTER NAME OF FILE rotation-trans op (.o) TO BE INVERTED:' 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 50 j=1,4 C50 read(2,1151)(M1(I,J),I=1,3) 50 read(2,*)(M1(I,J),I=1,3) 1151 FORMAT(3f14.7) CLOSE(UNIT=2) MM=3 NN=4 CALL MATPRN (M1,MM,NN) do 60 I=1,3 60 M1(i,i)=M1(i,i)-1 CALL MATPRN (M1,MM,NN) c invert matrix (top 3 rows) of operator N=3 CALL MATINV (M1,N,DET) write(6,*)'Det=',DET c ***note that since transpose of rotaion matrix is its inverse, we could skip matinv!****** c *** for now use real matinv for generality. c multiply translation by inverse matrix, take with negative sign as new translation do 100 i=1,3 m2(i,4)=0 do 100 j=1,3 100 m2(i,4)=m2(i,4)+m1(i,j)*m1(j,4) do 110 i=1,3 110 m1(i,4)=-m2(i,4) WRITE(6,*)' NCS OPERATOR:' WRITE(6,*)' ROTATE MATRIX' DO 1100 I=1,3 1100 WRITE(6,9) (M1(I,J),J=1,3) WRITE(6,*)' TRANSLATE ' WRITE(6,9) (M1(I,4),I=1,3) 9 FORMAT(' ',3F14.6) write(6,*)' this operator will be saved as an O datablock RT.O' WRITE(6,*)' Note this operator is for orthogpnal coord, A translation' WRITE(6,*)' You will have to convert for your cell. ' WRITE(6,*) 'ENTER NAME FOR output of resulting rot-trans op (.o):' READ(5,1200) filename open (unit=3,name=FILENAME,status='unknown') WRITE(3,*) '.SPACE_GROUP_operators r 12 (3f14.7)' do 1150 j=1,4 1150 write(3,1151)(M1(I,J),I=1,3) 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