CHARACTER*40 INFILE,FILENAME character*3 type, ID 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),y(3) character*40 string INTEGER*4 MM,NN,I,J c read operator into m1 WRITE(6,*) 'Find point on axis of a pure rotation RT-op.' WRITE(6,*) 'For arb RT-op, find point on axis of the pure' WRITE(6,*) 'rotation that could be made by adjusting z-trans' WRITE(6,*) 'to give no translation along axis' WRITE(6,*) 'This may not perform well if rot axis is near' WRITE(6,*) 'the x-y plane. see /text/rotorigin.txt for' WRITE(6,*) 'ideas for improvement.' 5 WRITE(6,*) 'ENTER NAME OF FILE rotation-trans op' 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) y(1)=0 c T2 - T1m23/m13 c ---------------------- = (y2) ****(2)**** c (m22- 1 - m12m23/m13) y(2)= -M1(2,4) + (M1(1,4)*M1(2,3))/M1(1,3) y(2)= y(2)/(M1(2,2)-(1+M1(1,2)*M1(2,3)/M1(1,3))) c y3= (T1 - m12y2)/m13 ****(1)**** y(3)= (-M1(1,4) - M1(1,2)*y(2))/M1(1,3) write(6,*) 'Y=',(y(i),i=1,3) c----------------------------------------------------------------- c---Now calculate vector along rotation axis of matrix: c subtract Identity matrix from rotation matrix c if original matrix took point to itself, this takes to zero c for point to go to zero, nust be orthoganl to the three rows of matrix c for this to have solution other than 0,0,0; the 3 rows coplanar. do 70 i=1,3 70 M1(I,I)=M1(I,I)-1.0 i=1 j=2 s=(m1(j,3)/m1(j,1)-m1(i,3)/m1(i,1))/(m1(j,2)/m1(j,1)-m1(i,2)/m1(i,1)) r=(-m1(i,2)/m1(i,1))*s + (-m1(i,3)/m1(i,1)) write(6,*) 'r, s, 1=',r, s, 1. write(6,*) 'Any multiple of this vector added to Y above' write(6,*) 'gives a point on the rotation axis.' iatom=101 type='C' ID='ROT' occ=0.0 b=0.0 write(6,*)'Points on rotation axis:' write(6,*)y(1)+scale*r,y(2)+scale*s,y(3)+scale*1. write(6,*)y(1)-scale*r,y(2)-scale*s,y(3)-scale*1. WRITE(6,51) iatom,type,ID,iatom,(y(I),i=1,3),occ,b do 140 i=10,1,-1 iatom=iatom+1 scale=float(i)/10. 140 WRITE(6,51) iatom,type,ID,iatom, . y(1)-scale*r,y(2)-scale*s,y(3)-scale*1.,occ,b do 150 i=1,10 iatom=iatom+1 scale=float(i)/10. 150 WRITE(6,51) iatom,type,ID,iatom, . y(1)+scale*r,y(2)+scale*s,y(3)+scale*1.,occ,b 51 format ('ATOM ',I5,2x,A1,3x,A3,' X',i4,4x,3f8.3,2f6.2) 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 cATOM 101xxxxxx A 101 0.000 -34.018 -7.215 0.00 0.00 cATOM 3103 O1D HEM A 2 15.332 -8.749 60.863 1.00 20.00 cATOM 103 C ROT A 103 15.875 -32.558 -8.215 0.00 0.00