CHARACTER*40 INFILE,FILENAME LOGICAL EFLG real*8 M1(3),M2(3,4),M3(3),M(3,4),ACC(3,4) REAL*8 P1(3),P2(3),T(3) character*40 string INTEGER*4 MM,NN,I,J 1200 FORMAT(A) c if crit is smaller than screw component of operator, won't converge but c keeps moveing by this vector along axis. Use this to correct translatio c component of operator write(6,5)'enter criterion for convergence (movement on operation):' read (5,*)crit c read operator into m2 6 WRITE(6,*) 'ENTER NAME OF rotation-trans op (.o):' READ(5,1200) INFILE IF (INFILE(:5).EQ.' ') STOP INQUIRE(FILE=INFILE, NAME=FILENAME,EXIST=EFLG) IF (.NOT.EFLG) GOTO 6 OPEN (UNIT=2,FILE=INFILE,STATUS='OLD',READONLY) c READ(2,*) '.SPACE_GROUP_operators r 12 (3f12.7)' READ(2,'(A40)') string do 60 j=1,4 60 read(2,*)(M2(I,J),I=1,3) CLOSE(UNIT=2) write(6,*) MM=3 NN=4 c CALL MATPRN (M2,MM,NN) IATOM=1 c read vector into m1 5 WRITE(6,*)'ENTER coordinates of point to be operated on:' 50 read(5,*,end=1160) (m1(i),I=1,3) c write(6,*) (m1(i),I=1,3) c CALL MATPRN (M1,MM,NN) c Multiply m1 matrix and transl by m2 matrix (first 3 cols), put im M3 55 write(6,9) 1,(m1(i),I=1,3) do 100 i=1,3 do 100 j=1,1 m3(i)=0 do 100 k=1,3 m3(i)=m3(i)+m2(i,k)*m1(k) C WRITE(6,*)I,J,K,M2(I,K),M1(K),M3(I) 100 CONTINUE c now add the translation of the second operator: do 110 i=1,3 110 m3(i)=m3(i)+m2(i,4) c WRITE(6,9) iatom,iatom,(M3(I),i=1,3),1.0,0.0 51 format (' ATOM ',I5,' N ALA A',i5,' ',3f8.3,2f6.2) WRITE(6,9) 2,(M3(I),i=1,3) 9 FORMAT(' 'i2,3F14.6) IATOM=IATOM+1 c now take midpoint between input point and output point as estimate of point on axis; repeat r=0. do 120 i=1,3 x=m3(i)-m1(i) r=r+x*x 120 m1(i)=m1(i)+x/2. r=sqrt(r) write(6,*)r if (r.gt.crit) goto 55 write(6,*) 1150 goto 50 1160 continue 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