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) character*40 fname pi=4*atan(1.0) c alpha helix has 100 degree rotation per radius, c 1.5 translation along helix, and c N atoms are 1.5 A? from axis dtheta = 100.*pi/180. dz=1.5 rh=1.5 C write(6,*)'What should we use for distance of NN nitrogen from axis?' C read(5,*)rh do 100 i=1,30 j=i-1 theta=j*dtheta b(i,1)=rh*cos(theta) b(i,2)=rh*sin(theta) b(i,3)=j*dz c write(6,401)(b(i,j),j=1,3) 100 continue c A has the coordinates of nitrogens of helix at arbitrary orientation c B has coordinates of nitrogens of helix along z axis, 1st at 1.5,0,0 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,*)' Enter name of file w coordinates of N to be fit:' read(5,'(A)') fname 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) c write(6,401) (a(i,j),j=1,3) 150 continue 151 format(30x,3f8.3) c50 format (4f10.3,3I5,f8.2,A15) RDI C50 format ('HETATM',I5,' N ALA A 1 ',3f8.3,2f6.2) c translate to bring MIDDLE atom Mid of both sets to 0,0,0: Mid=N/2 do 160 j=1,3 160 v(j)=a(Mid,j) do 170 i=1,n do 170 j=1,3 170 a(i,j)=a(i,j)-v(j) do 180 j=1,3 180 v(j)=b(Mid,j) do 190 i=1,n do 190 j=1,3 190 b(i,j)=b(i,j)-v(j) 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 (PDB helix)------- ------B (calc helix)----- ', 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 DO 640 I=1,N DO 640 J=1,3 640 A(I,J)=A(I,J)-VS(J) GOTO 611 650 CONTINUE 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