real*4 A(30,3), B(30,3),ATA(3,3),ATB(3,3),R(3,3),v(3),VS(3),temp(3,3) real*4 M(3,3),SQM(3,3),coma(3),comb(3) character*40 fname 100 continue c A has the coordinates of set A (maybe in pixels not A) c B has the coordinates of set B (maybe in pixels not A) c ATA is ATranspose * A, or its inverse c ATB is Atranspose * B c R is rotation matrix we are solving for, which will take A to B write(6,*)' Program to find rot-trans operator to take points' write(6,*)' in set A to points in set B.' write(6,*)' each set in PDB file, atom recs only, n as first line' write(6,*) write(6,*)' Enter name of PDB file w coordinates of A to be fit:' read(5,'(A)') fname write(6,*)'A:' open (unit=2,name=fname,status='old') read(2,*) n IF (N.GT.30) N=30 do 150 i=1,n read(2,151) (a(i,j),j=1,3) write(6,401) (a(i,j),j=1,3) 150 continue 151 format(30x,3f8.3) close (unit=2) write(6,*)' Enter name of file w coordinates of B:' read(5,'(A)') fname open (unit=2,name=fname,status='old') write(6,*)'B:' read(2,*) j IF (j.NE.n) write(6,*)'using',n,' records not',j do 160 i=1,n read(2,151) (b(i,j),j=1,3) write(6,401) (b(i,j),j=1,3) 160 continue close (unit=2) write(6,*)'Enter matrix to bring these coord to', . ' orthogonal A coordinates:' do 170 i=1,3 170 read(5,*)(ata(i,j),j=1,3) do 190 i=1,n do 180 j=1,3 v(j)=0 do 180 k=1,3 180 v(j)=v(j)+ata(j,k)*a(i,k) do 190 j=1,3 190 a(i,j)=v(j) do 1210 i=1,n do 1200 j=1,3 v(j)=0 do 1200 k=1,3 1200 v(j)=v(j)+ata(j,k)*b(i,k) do 1210 j=1,3 1210 b(i,j)=v(j) write(6,*)'A after orthog:' do 1211 i=1,N 1211 write(6,401) (a(i,j),j=1,3) write(6,*) write(6,*)'B after orthog:' do 1212 i=1,N 1212 write(6,401) (b(i,j),j=1,3) write(6,*) c check rms difference: C ZERO VS(),ssx2,TE FOR ACCUMULATING ERROR 1611 DO 1612 I=1,3 1612 VS(I)=0 ssx2=0 te=0 write(6,613) do 1620 i=1,n sx2=0 do 1615 j=1,3 v(j)=a(i,j)-b(i,j) VS(J)=VS(J)+V(J) 1615 sx2=sx2+v(j)*v(j) ssx2=ssx2+sx2 sx2=sqrt(sx2) 1620 write (6,621) (A(i,j),j=1,3),(B(i,j),j=1,3),(v(j),j=1,3),sx2 write(6,*)' rh=',rh,' rms residual:',sqrt(ssx2/N) DO 1630 J=1,3 vs(j)=vs(j)/n 1630 TE=TE+VS(J)*vs(j) TE=SQRT(TE) WRITE(6,*)' TRANSLATIONAL ERROR, A-B:',(VS(J),J=1,3) write(6,*)' SCALAR:',TE,' A.' c get center of mass of a: do 1250 j=1,3 coma(j)=0 do 1250 i=1,n 1250 coma(j)=coma(j)+a(i,j) c get center of mass of b: do 1260 j=1,3 comb(j)=0 do 1260 i=1,n 1260 comb(j)=comb(j)+b(i,j) do 1270 j=1,3 coma(j)=coma(j)/n 1270 comb(j)=comb(j)/n c subtract com from each: do 1280 i=1,n do 1280 j=1,3 a(i,j)=a(i,j)-coma(j) 1280 b(i,j)=b(i,j)-comb(j) write(6,*)'A after orthog, centering :' do 1281 i=1,N 1281 write(6,401) (a(i,j),j=1,3) write(6,*) write(6,*)'B after orthog, centering :' do 1282 i=1,N 1282 write(6,401) (b(i,j),j=1,3) write(6,*) c form ATA, a 3x3 matrix of the dp of columns of A do 200 i=1,3 do 200 k=1,3 ata(i,k)=0 do 200 j=1,n 200 ata(i,k)=ata(i,k)+a(j,i)*a(j,k) !+At(i,j)*A(j,k) write(6,*)'ATA:' call matprn(ata,3,3) write(6,*) c form ATB, product of A-transpose and B do 250 i=1,3 do 250 k=1,3 atb(i,k)=0 do 250 j=1,n c write(6,*) i,k,j,a(j,i),b(j,k),atb(i,k) 250 atb(i,k)=atb(i,k)+a(j,i)*b(j,k) !+At(i,j)*A(j,k) write(6,*)'ATB:' call matprn(atb,3,3) write(6,*) c save in temp and invert ATA: c temp=ATA do 252 i=1,3 do 252,j=1,3 252 temp(i,j)=ATA(i,j) call matinv(ATA,3,det) write(6,*)'ATA-inv:' call matprn(ata,3,3) write(6,*) c test inversion and least squares inversion: ATA-inv * ATA =1 write(6,*)' ATA-inv * ATA:' do 260, i=1,3 do 255, k=1,3 v(k)=0 do 255 j=1,3 255 v(k)=v(k)+ATA(i,j)*temp(j,k) 260 write(6,401) (v(k),k=1,3) c R = ATA-inv*ATB do 300 i=1,3 do 300 k=1,3 R(i,k)=0 do 300 j=1,3 300 R(i,k)=R(i,k)+ATA(i,j)*ATB(j,k) c print R matrix on screen: write(6,*) write(6,*)'R, = ATA-inv * ATB:' CALL MATPRN(R,3,3) 401 format(' ',2(3f15.8,3x)) C DECOMPOSE GENERAL MATRIX R INTO ORTHONORMAL R AND DIAGONAL SQM CALL DECOMP(R,SQM,3) c print R matrix on screen: write(6,*) write(6,*)'R, ORTHONORMAL COMPONENT OF PREVIOUS:' CALL MATPRN(R,3,3) c write R matrix to file write(6,*)'Writing TRANSPOSE of R to file rmatrix.mat' write(6,*)' i.e. matrix to pre-multiply pdb to' write(6,*)' align helix along +z.' write(6,*)' (R is to post-multiply PDB)' open(unit=3,name='rmatrix.mat',status='new') do 500 i=1,3 c500 write(3,401) (R(i,j),j=1,3) 500 write(3,401) (R(j,i),j=1,3) close(unit=3) c apply rot matrix to A and calc residuals A2 = A * R (mult by R on right!) do 610 i=1,n do 600 j=1,3 v(j)=0 do 600 k=1,3 600 v(j)=v(j)+A(i,k)*R(k,j) do 610 j=1,3 610 A(i,j)=v(j) C ZERO VS(),ssx2,TE FOR ACCUMULATING ERROR 611 DO 612 I=1,3 612 VS(I)=0 ssx2=0 te=0 write(6,613) 613 format('------A------------------- ------B------------------ ', 6 '----------diff---------- --Angst--') do 620 i=1,n sx2=0 do 615 j=1,3 v(j)=a(i,j)-b(i,j) VS(J)=VS(J)+V(J) 615 sx2=sx2+v(j)*v(j) ssx2=ssx2+sx2 sx2=sqrt(sx2) 620 write (6,621) (A(i,j),j=1,3),(B(i,j),j=1,3),(v(j),j=1,3),sx2 621 format(' ',4(3f8.3,3x)) write(6,*)' rh=',rh,' rms residual:',sqrt(ssx2/N) DO 630 J=1,3 vs(j)=vs(j)/n 630 TE=TE+VS(J)*vs(j) TE=SQRT(TE) WRITE(6,*)' TRANSLATIONAL ERROR, A-B:',(VS(J),J=1,3) write(6,*)' SCALAR:',TE,' A.' C IF TE GE .1 A, SUBTRACT IT FROM ALL ATOMS AND RECALCULATE ERRORS IF (TE.LT.0.1) GOTO 650 write(6,*) ' an additional offset is being subtracted. This is NOT' write(6,*)' included in the output operator!' DO 640 I=1,N DO 640 J=1,3 640 A(I,J)=A(I,J)-VS(J) GOTO 611 650 CONTINUE c get final offset= comb+[M](coma) do 670 i=1,3 v(i)=0 do 670 j=1,3 670 v(i)=v(i)+R(j,i)*coma(j) WRITE(6,*)' RT OPERATOR:' WRITE(6,*)' ROTATE MATRIX' DO 1100 I=1,3 1100 WRITE(6,9) (R(J,I),J=1,3) !R IS TRANSPOSE OF STANDARD ROTATION MATRIX! WRITE(6,*)' TRANSLATE ' WRITE(6,9) (COMB(I)-V(I),I=1,3) !COMB - [M]*COMA 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,'(A)') fname c write(6,*)'xx',filename,'xx' open (unit=3,name=FNAME,status='unknown') WRITE(3,*) '.SPACE_GROUP_operators r 12 (3f12.7)' do 1150 j=1,3 1150 write(3,1151)(R(J,I),I=1,3) WRITE(3,1151) (COMB(I)-V(I),I=1,3) !COMB - [M]*COMA 1151 FORMAT(3f12.7) END c----------------------------------------------------------------- SUBROUTINE DECOMP(R,S1,N) C DECOMPOSE GENERAL MATRIX R INTO ORTHONORMAL R AND DIAGONAL S2 REAL*4 R(3,3),M(3,3),S1(3,3),S2(3,3),S1I(3,3),V(3) c test orthonormality of "rotation" matrix: write(6,*) write(6,*)'Rt*R, test orthonormality:' do 510 i=1,3 do 505 k=1,3 M(I,K)=0 do 505 j=1,3 505 M(I,k)=M(I,k)+R(j,i)*R(j,k) !RT(i,j)*R(j,k) 510 write(6,401) (M(I,k),k=1,3) 401 format(' ',2(3f15.8,3x)) C FIND S, THE SYMMETRIC? SQRT OF M: SQM*SQM=M, SQM=SQM-INV * M C ITERATIVE METHOD, S1 IS APPROX, S2 WILL BE BETTER APPROX C FIRST APPROX: S1 IS DIAGONAL, WITH ELEMENTS = SQRT(DIAG OF M) DO 100 I=1,3 DO 90 J=1,3 90 S1(I,J)=0 100 S1(I,I)=SQRT(M(I,I)) 110 CONTINUE C COPY S1 TO S1I (= S1-INV) AND INVERT: DO 120 I=1,3 DO 120 J=1,3 120 S1I(I,J)=S1(I,J) CALL MATINV(S1I,3,DET) C MULTIPLY M BY S1I AND PUT RESULT IN S2 DO 140 I=1,3 DO 140 K=1,3 S2(I,K)=0 DO 140 J=1,3 140 S2(I,K)=S2(I,K)+S1I(I,J)*M(J,K) C MEASURE DIFFERENCE BETWEEN S1 AND S2, AND PUT AVERAGE INTO S2 SSQE=0 DO 160 I=1,3 DO 160 J=1,3 X=S1(I,J)-S2(I,J) SSQE=SSQE+X*X 160 S2(I,J)=S2(I,J)+X/2 WRITE(6,*)'SSQE=',SSQE C IF ERROR OK, GOT TO END AND RETURN. OTHERWISE COPY S2 TO S1 AND REPEAT. IF (SSQE.LT.1E-12) GOTO 900 DO 200 I=1,3 DO 200 J=1,3 200 S1(I,J)=S2(I,J) GOTO 110 C MULTIPLY R BY S-INV ON THE LEFT, PUT RESULT IN R (temporarily in S2) 900 DO 300 I=1,3 DO 300 K=1,3 S2(I,K)=0 DO 300 J=1,3 300 S2(I,K)=S2(I,K)+R(I,J)*S1I(J,K) C TEST IF R*S = ORIG MATRIX (R) SSQE=0 DO 410 I=1,3 DO 400 K=1,3 S1I(I,K)=0 DO 400 J=1,3 400 S1I(I,K)=S1I(I,K)+S2(I,J)*S1(J,K) DO 405 K=1,3 V(K)=R(I,K)-S1I(I,K) 405 SSQE=SSQE+V(K)*V(K) 410 WRITE(6,621)(R(I,K),K=1,3),(S1I(I,K),K=1,3),(V(K),K=1,3) 621 format(' ',4(3f9.4,3x)) WRITE(6,*)' RMSERROR BETWEEN ORIG MAT AND S*R =',SQRT(SSQE) C COPY R-MATRIX FROM S2 TO R. S MATRIX STILL IN S1 (BETTER WAS IN S2, BUT R CALC USING INV OF S1) DO 500 I=1,3 DO 500 J=1,3 500 R(I,J)=S2(I,J) c test orthonormality of "rotation" matrix: write(6,*) write(6,*)'Rt*R, test orthonormality:' do 610 i=1,3 do 605 k=1,3 M(I,K)=0 do 605 j=1,3 605 M(I,k)=M(I,k)+R(j,i)*R(j,k) !RT(i,j)*R(j,k) 610 write(6,401) (M(I,k),k=1,3) RETURN END C------------------------------------------------------------------ SUBROUTINE MATINV (ARRAY, NORDER, DET) real*4 ARRAY, AMAX, SAVE 2 DIMENSION ARRAY(3,3), 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 (ABS (AMAX) - ABS (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 SUBROUTINE MATPRN(ARRAY,M,N) C DOUBLE PRECISION ARRAY REAL*4 ARRAY(3,3) DO 10 I=1,M 10 WRITE (6,*) (ARRAY(I,J),J=1,N) RETURN END