C PROTEIN DATA BANK SOURCE CODE THEOD C AUTHOR L.LEBIODA C ENTRY DATE. 1/82 UNSUPPORTED C LAST REVISION. 1/82 C PURPOSE. MEASURE COORDINATES WITH C PURPOSE. THEODOLITE. C LANGUAGE. FORTRAN IV C C C PROGRAM THEOD(INPUT,OUTPUT,TAPE1,TAPE2,TAPE5=INPUT,TAPE6=OUTPUT) C C THIS IS AN INTERACTIVE PROGRAM TO GUIDE THE MEASUREMENT OF THE C COORDINATES OF A KENDREW MODEL WITH TWO THEODOLITES (SURVEYORS C TRANSITS) AND TO CALCULATE CARTESIAN COORDINATES FROM THE MEASU- C REMENTS. C C THEODOLITES USUALLY CAN BE RENTED FROM DEALERS OF SURVEYING INSTRU- C MENTS. INSTRUMENTS WITH AN ACCURACY OF 1 MINUTE OR 0.01 GRAD CAN BE C USED; CONVENIENCE OF ANGLE READING IS VERY IMPORTANT AND THEREFORE C TRANSITS ARE NOT RECOMENDED. THE USE OF THEODOLITES IS SIMPLE AND C NO TRAINING IS NECESSARY EXCEPT FOR SOME INSTRUCTIONS FROM THE C DEALER. C THE PRESENT VERSION IS FOR INPUT IN GRADS (PI=200 GRADS), CARDS FOR C INPUT IN DEGREES AND MINUTES HAVE BEEN COMMENTED OUT AND CAN BE C REACTIVATED (IN THEOD, MEASURE AND RADIAN). C AZIMUTHAL ANGLE (FI) SHOULD BE READ COUNTERCLOCKWISE, OTHERWISE C THEOD AND COOR SHOULD BE CHANGED AS COMMENTED. C WHEN STARTING THE MEASUREMENT, THE THEODOLITES SHOULD BE POSITIONED C RELATIVELY CLOSE TO THE MEASURED MODEL BUT SO AS TO BE ABLE TO C FOCUS THE SCOPE ON THE CLOSEST OBJECTS TO BE MEASURED. THE DISTANCE C BETWEEN THE THEODOLITES SHOLD ALSO BE GREATER THAN THE MINIMUM C FOCAL LENGTH. A GOOD SETUP IS ONE CLOSE TO A RIGHT ISOSCALES C TRIANGLE WITH THE RIGHT ANGLE AT THE MODEL. C TO SCALE THE SETUP, TWO POINTS WITH PRECISELY KNOWN DISTANCE, C ABOUT 100 CM, SHOULD BE AVAILABLE (ENDS OF MACHINISTS SCALE). C C THE MEASUREMENT CONSISTS OF AIMING THE CROSSHAIRS OF BOTH TELE- C SCOPES AT THE OBJECT (ESTIMATED POSITON OF AN ATOM) AND INPUTING C THE ANGLES AND THE OBJECT IDENTIFIER. IT IS USEFUL, ESPECIALLY C IN THE CROWDED PARTS OF THE MODEL, TO POINT THE MEASURED ATOM C WITH A PIECE OF COLORED MASKING TAPE. C C IT IS CONVENIENT TO INTRODUCE A COORDINATE SYSTEM BASED ON THE C ELECTRON DENSITY MAP RATHER THAN ON THE THEODOLITE SETUP. C TO DO THIS, A FEW STANDARDS (AT LEAST 5 BUT NO MORE THAN 30) C SHOULD BE MOUNTED IN THE RICHARDS BOX SO THAT THEIR MAP COORDINA- C TES ARE EASYLY DETERMINED (THEY ARE REFERRED TO IN THE PROGRAM AS C TARGET COORDINATES). THESE STANDARDS ARE SUBSEQUENTLY MEASURED C WITH THE THEODOLITES AND ISOMETRIC TRANSFORMATION TO THE MAP C COORDINATE SYSTEM IS CALCULATED BY THE LEAST SQURES METHOD. C WITH THIS PROCEDURE IT IS ALSO POSSIBLE TO MERGE MEASUREMENTS C DONE WITH DIFFERENT THEODOLITE SETUPS. DIFFERENT ARRANGEMENTS C OF THEODOLITES CAN RESULT FROM ACCIDENTAL PUSHING OF ONE OF C THE INSTRUMENTS OR BE INTRODUCED ON PURPOSE TO AVOID AN OBSCURED C VIEW OF A PART OF THE MODEL. THIS, HOWEVER, IS SELDOM NECEESSARY. C ALTHOUGH AN ATOM IS SOMETIMES OBSCURED BY INTERVENING MODEL C PARTS, IN MOST CASES IT IS STILL POSSIBLE TO AIM CROSSHAIRS OF C THE TELESCOPE BY NOTING TWO BONDS THAT INTERSECT AT THE DESIRED C ATOMIC POSITION AND VISUALLY EXTRAPOLATING THEIR (NON-VISIBLE) C INTERSECTION POINT. THE RELIBILITY OF SUCH A MEASUREMENT IS VERFIED C BY AN UNCERTAINITY PARAMETER CALCULATED BY THE PROGRAM. IF AIMING C WAS POOR IT CAN BE REPEATED. C C FINAL COORDINATES ARE WRITTEN ON TAPE1. TAPE2 IS USED FOR INTER- C MEDIATE RESULTS, INFORMATION ABOUT COORDINATE SYSTEMS AND C THE TRANSFORMATION USED. IF THE THEODOLITES HAVE NOT BEEN MOVED C THIS INFORMATION CAN BE USED TO RESTART MEASUREMENT WITHOUT C REMEASURING THE STANDARDS DETERMINING THE COORDINATE SYSTEM C (UNDISTURBED RESTART). C C IF THE ATOM IDENTIFIER IS RESIDUE NAME, RESIDUE NUMBER AND ATOM NAME C AS FOR EXAMPLE THR129CD2 , THEN THE OUTPUT FILE IS IN THE PROTEIN C DATA BANK FORMAT. C TO EXIT FROM THE PROGRAM INPUT "END" AS ATOM IDENTIFIER. C C PROGRAM IS IN VERY BASIC FORTRAN AND SHOULD BE MACHINE INDEPENDENT C EXCEPT FROM THE CARD 'PROGRAM'. C C C THIS PROGRAM WAS WRITTEN BY C C LUKASZ LEBIODA C DEPARTMENT OF CHEMISTRY, JAGIELLONIAN UNIVERSITY, KRAKOW, POLAND C WHILE ON LEAVE OF ABSENCE AT THE DEPARTMENT OF CHEMISTRY, MICHIGAN C STATE UNIVERSITY, EAST LANSING. C C SOME SUBROUTINES FROM PROGRAM HOMOMGR BY M.G.ROSSMANN WERE ADOPTED C TO THIS PROGRAM. C C THIS WORK WAS SUPPORTED BY GRANT GM 26360 FROM N.I.H. C TO C ALEXANDER TULINSKY C DEPARTMENT OF CHEMISTRY, MICHIGAN STATE UNIVERSITY, EAST LANSING. C C C DIMENSION A(3),X(3),XX(3),XL(3),XP(3),XT(3,30),XM(3,30),XR(3,30), +R(3,4),Y(3),THE(3),FI(2),CHI(2) C THE FOLLOWING CARD SHOULD BE REACTIVATED FOR INPUT IN DEGREES C AND MINUTES C INTEGER TFI1(2),TCHI1(2),TFI2(2) COMMON/L/ DL,FI1T,FI2T,TAU,NEG LOGICAL NEG DATA R/1.0,3*0.0,1.0,3*0.0,1.0,3*0.0/ DATA Y/3*0.0/ DATA PI/3.14159266/ WRITE(6,100) WRITE(6,101) 100 FORMAT(33H WELCOME TO PROGRAM THEOD +,/,48H ALLIGN THEODOLITES IN THEIR HORIZONTAL PLANES. +,/,44H WHEN STANDING BETWEEN THEM FACING THE MODEL +,/,44H THE ONE AT LEFT IS REFERRED AS NUMBER ONE ,/) C THE LAST CARD SHOULD BE COMMENTED OUT AND SUBSTITUTED WITH THE C NEXT FOR CLOCKWISE READOUT OF THEODOLITES C +,/,45H THE ONE AT RIGHT IS REFFERED AS NUMBER ONE ,/) 101 FORMAT(36H POINT THEODOLITES EACH AT THE OTHER +,/,46H INPUT FI1, CHI1 (FROM THE FIRST THEOD.), FI2 ,/) READ (5,*) TFI1,TCHI1,TFI2 11 FORMAT(4F10.4) FI1T=RADIAN(TFI1) FI2T=RADIAN(TFI2) TAU=RADIAN(TCHI1)-1.57079633 NEG=.FALSE. IF(TAU) 41,42,42 41 NEG=.TRUE. TAU=-TAU 42 CONTINUE WRITE(6,50) 50 FORMAT(33H DO YOU WANT UNDISTURBED RESTART? ,/, +11H YES OR NO ,/) READ(5,1)Q IF(Q.NE.1HY)GO TO 90 51 WRITE(6,52) 52 FORMAT(36H INPUT DISTANCE BETWEEN THEODOLITES ,/) READ(5,*) DL WRITE(6,53) 53 FORMAT(34H INPUT ROTATION MATRIX AS PRINTED ,/) READ(5,*)((R(I,J),J=1,4),I=1,3) WRITE(6,105)DL WRITE(6,54) 54 FORMAT(20H ROTATION MATRIX IS ,/) WRITE(6,11)((R(I,J),J=1,4),I=1,3) WRITE(6,55) 55 FORMAT(25H IS IT O.K.? (YES OR NO) ,/) READ(5,1)Q IF(Q.NE.1HY) GO TO 51 WRITE(2,105)DL WRITE(2,54) WRITE(2,11) ((R(I,J),J=1,4),I=1,3) WRITE(2,110) GO TO 30 90 WRITE(6,102) 102 FORMAT(42H INPUT VERY CRUDE ESTIMATE OF THE DISTANCE +,/,37H BETWEEN THEODOLITES IN CENTIMETERS ,/) READ(5,*)DL DL=0.5*DL WRITE(6,103) 103 FORMAT(51H MEASURE TWO POINT WITH KNOWN DISTANCE BETWEEN THEM ,/) CALL MEASURE (XL,A,END,XL) CALL MEASURE (XP,A,END,XL) P=0.0 DO 7 I=1,3 Q=XL(I)-XP(I) 7 P=P+Q*Q P=SQRT(P) WRITE(6,104) 104 FORMAT(39H INPUT THE DISTANCE BETWEEN THE POINTS , +26HIN THE RICHARDS BOX SCALE ,/) READ(5,*) D DL=DL*D/P WRITE(2,109) WRITE(6,109) 109 FORMAT(38H FROM NOW DISTANCES ARE IN ANGSTREMS ,/) WRITE(2,105) DL WRITE(6,105) DL 105 FORMAT(1X,32H DISTANCE BETWEEN THEODOLITES IS ,F10.4,/) 12 WRITE(6,106) 106 FORMAT(40H IF YOU WANT TO SET THE REFERENCE SYSTEM +,/,54H BASED ON KNOWN COORDINATES OF N POINTS INPUT N.GT.4 ,/) WRITE(2,110) 110 FORMAT(40H0 ATOM FI1 CHI1 FI2 +,56H CHI2 X Y Z UNCERT DISTANCE ) READ(5,*) NP IF(NP.LE.4) GO TO 30 DO 20 N=1,NP WRITE(6,107) 107 FORMAT(24H INPUT KNOWN COORDINATES ,/) READ (5,*) (XT(I,N),I=1,3) CALL MEASURE(X,A,END,Y) IF(END.EQ.3HEND) GO TO 99 DO 21 M=1,3 Y(M)=X(M) 21 XM(M,N)=X(M) 20 CONTINUE CALL LINEAR(XT,XM,R,THE,NP) CALL REFEUL(XT,XM,R,THE,NP) CALL ROTATE(XM,XR,R,NP) CALL COMPARE(XT,XR,NP) WRITE(6,108) 108 FORMAT(36H DO YOU ACCEPT THIS TRANSFORMATION +,/,12H YES OR NO ,/) READ(5,1)Q 1 FORMAT(A1) IF(Q.NE.1HY) GO TO 12 30 CONTINUE CALL MEASURE(X,A,END,Y) IF(END.EQ.3HEND) GO TO 99 DO 31 I=1,3 Y(I)=X(I) XX(I)=R(I,4) DO 32 J=1,3 XX(I)=XX(I)+R(I,J)*X(J) 32 CONTINUE 31 CONTINUE C C THIS WRITE STATEMENT IS FOR OUTPUT FILE C WRITE(1,300) A(3),A(2),A(1),XX 300 FORMAT(4HATOM,9X,A3,1X,A3,3X,A3,4X,3F8.3) GO TO 30 99 ENDFILE 1 ENDFILE 2 STOP END SUBROUTINE MEASURE(X,A,END,Y) DIMENSION A(3),X(3),Y(3),FI(2),CHI(2) C C THE FOLLOWING CARD IS FOR INPUT IN DEGREES AND MINUTES C INTEGER FI1(2),CHI1(2),FI2(2),CHI2(2) COMMON/L/ DL,FI1T,FI2T,TAU,NEG 100 WRITE(6,10) 10 FORMAT(37H POINT THE THEODOLITES AT THE OBJECT +,/,34H INPUT POINT IDENTIFICATION (3A3) ,/) READ(5,1) A 1 FORMAT(3A3) END=A(1) IF(END.EQ.3HEND)RETURN WRITE(6,11) 11 FORMAT (48H INPUT FI1, CHI1, FI2, CHI2 ,/) READ (5,*) FI1,CHI1,FI2,CHI2 12 FORMAT (4F10.4) FI(1)=RADIAN(FI1) FI(2)=RADIAN(FI2) CHI(1)=RADIAN(CHI1) CHI(2)=RADIAN(CHI2) CALL COOR(FI,CHI,X,S) WRITE(6,2) A,FI1,CHI1,FI2,CHI2 C C THE FOLLOWING CARD FOR INPUT IN DEGREES AND MINUTES C 2 FORMAT(1X,3A3,4(1X,I5,I3,1X),3F8.3,3X,3F8.3) C THE FOLLOWING CARD SHOULD BE COMMENTED OUT FOR INPUT IN DEGREES 2 FORMAT(1X,3A3,4F10.4,3F8.3,3X,3F8.3) WRITE(6,7) X 7 FORMAT(6H XYZ=,3F8.3) WRITE(6,3) 0.5*S 3 FORMAT(15H UNCERTAINITY=,F8.3,25H IS IT O.K.? (YES OR NO) ,/) 102 READ(5,4) Q 4 FORMAT(A1) IF(Q.EQ.1HN) GO TO 100 IF(Q.EQ.1HY) GO TO 101 WRITE(6,5) 5 FORMAT(13H YES OR NO? ,/) GO TO 102 101 CONTINUE D=0.0 DO 20 I=1,3 DD=X(I)-Y(I) 20 D=D + DD*DD D=SQRT(D) WRITE(6,6)D 6 FORMAT(38H DISTANCE FROM THE LAST MEASURED ATOM ,F8.3,/) WRITE(2,2) A,FI1,CHI1,FI2,CHI2,X,S,D RETURN END C C C SUBROUTINE COOR(FI,CHI,X,S) COMMON/L/ DL,FI1T,FI2T,TAU,NEG DIMENSION FI(2),CHI(2),A(2),B(2),C(2),X(3) DIMENSION P(3),Q(3) LOGICAL NEG DATA PI/3.1415926/ FI(1)=FI(1)-FI1T IF(FI(1).LT.0.0) FI(1)=FI(1) + 2.0*PI FI(2)=FI(2)-FI2T + PI IF(FI(2).GT.PI) FI(2)=FI(2) - 2.0*PI DO 1 I=1,2 AI=0.5*(CHI(I)-TAU) BI=0.5*(CHI(I)+TAU) SA=SIN(AI) SB=SIN(BI) CA=COS(AI) CB=COS(BI) IF(NEG) GO TO 8 CFI=COS(FI(I)) COTA=1.0/TAN(0.5*FI(I)) ATS=ATAN(COTA*SA/SB) ATC=ATAN(COTA*CA/CB) FI(I)=PI-ATS-ATC GO TO 9 8 CONTINUE C THIS FOR NEGATIVE TAU AL=PI-FI(I) CFI=COS(AL) COTA=1.0/TAN(0.5*AL) ATC=ATAN(COTA*CA/CB) FI(I)=ATAN(COTA*SA/SB) + ATC C 9 CHI(I)=ACOS(COS(TAU)*COS(CHI(I))+SIN(TAU)*SIN(CHI(I))*CFI) 1 CONTINUE C C HERE END OF POLAR COORDINATES TRANSFORMATION C DO 2 I=1,2 A(I)=COS(FI(I))*SIN(CHI(I)) B(I)=SIN(FI(I))*SIN(CHI(I)) C C THE FOLLOWING CARD FOR CLOCKWISE READING OF THE AZIMUTHAL ANGLE C B(I)=-B(I) C(I)=COS(CHI(I)) 2 CONTINUE AA= A(1)*A(2) + B(1)*B(2) + C(1)*C(2) AB= 1.0 - AA*AA U= DL*(A(1) -AA*A(2))/AB W= DL*(AA*A(1) -A(2))/AB P(1)=A(1)*U P(2)=B(1)*U P(3)=C(1)*U Q(1)=A(2)*W + DL Q(2)=B(2)*W Q(3)=C(2)*W S=0.0 DO 3 I=1,3 X(I)=(P(I) + Q(I))*0.5 S=S + (P(I) - Q(I))**2 3 CONTINUE S=SQRT(S) RETURN END C C C C THE FOLLOWING 5 CARDS FOR INPUT IN DEGREES AND MINUTES C FUNCTION RADIAN(DEGMIN) C INTEGER DEGMIN(2) C IF (DEGMIN(1)) 1,2,2 C 1 DEGMIN(2)=-DEGMIN(2) C 2 RADIAN=0.01745329*(FLOAT(DEGMIN(1))+0.016666667*FLOAT(DEGMIN(2))) C THE FOLLOWING 2 CARDS FOR INPUT IN GRADS FUNCTION RADIAN (GRADS) RADIAN=0.01570796*GRADS RETURN END C C C SUBROUTINE LINEAR (XT,XM,R,THE,NP) DIMENSION XT(3,30),XM(3,30),AA(4,4),VEC(4,3),R(3 1,4), THE(3), DIFF(3), CG(2,3), X1(3), X2(3) DO 101 I=1,2 DO 101 J=1,3 101 CG(I,J)=0.0 NCOUNT=0 DO 103 N=1,NP DO 102 J=1,3 CG(1,J)=CG(1,J)+XT(J,N) CG(2,J)=CG(2,J)+XM(J,N) 102 CONTINUE NCOUNT=NCOUNT+1 103 CONTINUE DO 104 I=1,2 DO 104 J=1,3 104 CG(I,J)=CG(I,J)/NCOUNT DO 105 I=1,4 DO 105 J=1,4 105 AA(I,J)=0.0 DO 106 I=1,4 DO 106 J=1,3 106 VEC(I,J)=0.0 DO 110 N=1,NP DO 107 I=1,3 X1(I)=XT(I,N) -CG(1,I) X2(I)=XM(I,N) -CG(2,I) 107 CONTINUE DO 109 I=1,3 DO 108 J=1,3 AA(I,J)=AA(I,J)+X2(I)*X2(J) VEC(I,J)=VEC(I,J)+X1(J)*X2(I) 108 CONTINUE AA(I,4)=AA(I,4)+X2(I) AA(4,I)=AA(I,4) VEC(4,I)=VEC(4,I)+X1(I) 109 CONTINUE AA(4,4)=AA(4,4)+1 110 CONTINUE CALL SIMINV (AA,4) DO 113 I=1,3 DO 112 J=1,4 R(I,J)=0.0 DO 111 K=1,4 R(I,J)=R(I,J)+AA(J,K)*VEC(K,I) 111 CONTINUE 112 CONTINUE 113 CONTINUE CALL ROTANG (R,THE) DO 114 I=1,3 114 DIFF(I)=0.0 NCOUNT=0 DO 117 N=1,NP DO 116 I=1,3 SUB=XT(I,N) DO 115 J=1,3 SUB=SUB-R(I,J)*XM(J,N) 115 CONTINUE DIFF(I)=DIFF(I)+SUB 116 CONTINUE NCOUNT=NCOUNT+1 117 CONTINUE DO 118 I=1,3 R(I,4)=DIFF(I)/NCOUNT 118 CONTINUE RETURN C END C C C SUBROUTINE ROTATE(XM,XR,R,NP) DIMENSION XM(3,30),XR(3,30),R(3,4) C C ROTATES MEASURED POINTS ONTO TARGET C WRITE (6,105) WRITE (2,105) WRITE (6,106) ((R(I,J),J=1,4),I=1,3) WRITE (2,106) ((R(I,J),J=1,4),I=1,3) DO 4 K=1,NP DO 2 I=1,3 XR(I,K)=R(I,4) DO 1 J=1,3 XR(I,K)=XR(I,K)+XM(J,K)*R(I,J) 1 CONTINUE 2 CONTINUE 4 CONTINUE RETURN 105 FORMAT(45H ROTATION MATRIX APPLIED TO MEASURED POINTS ,/) 106 FORMAT (2X,4F10.4) END C C C SUBROUTINE COMPARE(XT,XR,NP) DIMENSION XT(3,30), XR(3,30) WRITE(2,103) WRITE(6,103) DD2=0.0 DD=0.0 DO 2 N=1,NP D=0.0 DO 1 J=1,3 D=D + (XT(J,N) - XR(J,N))**2 1 CONTINUE DD2=DD2 + D D=SQRT(D) DD=DD + D WRITE(2,104) N,D,(XT(J,N),J=1,3),(XR(J,N),J=1,3) WRITE(6,104) N,D,(XT(J,N),J=1,3),(XR(J,N),J=1,3) 2 CONTINUE 103 FORMAT(17H POINT DISTANCE ,T30,6HTARGET,T55,8HMEASURED ,/) DD2=SQRT(DD2)/NP WRITE(2,105) DD2,DD/NP WRITE(6,105) DD2,DD/NP RETURN 104 FORMAT(2X,I5,F9.2,4X,3F8.2,4X,3F8.2) 105 FORMAT(20H DEVIATION STANDARD ,F6.2,5X,7HAVERAGE,F6.2) END SUBROUTINE CALC (THE,FN,DER,N) C C GIVEN EULERIAN ANGLES, CALCULATE COEFFICIENTS AND DERIVATIVES C DIMENSION THE(3), FN(9), DER(3,9) C1=COS(THE(1)) C2=COS(THE(2)) C3=COS(THE(3)) S1=SIN(THE(1)) S2=SIN(THE(2)) S3=SIN(THE(3)) FN(1)=-S1*C2*S3+C1*C3 FN(2)=C1*C2*S3+S1*C3 FN(3)=S2*S3 FN(4)=-S1*C2*C3-C1*S3 FN(5)=C1*C2*C3-S1*S3 FN(6)=S2*C3 FN(7)=S1*S2 FN(8)=-C1*S2 FN(9)=C2 IF (N.EQ.0) RETURN DER(1,1)=-C1*C2*S3-S1*C3 DER(2,1)=S1*S2*S3 DER(3,1)=-S1*C2*C3-C1*S3 DER(1,2)=-S1*C2*S3+C1*C3 DER(2,2)=-C1*S2*S3 DER(3,2)=C1*C2*C3-S1*S3 DER(1,3)=0. DER(2,3)=C2*S3 DER(3,3)=S2*C3 DER(1,4)=-C1*C2*C3+S1*S3 DER(2,4)=S1*S2*C3 DER(3,4)=S1*C2*S3-C1*C3 DER(1,5)=-S1*C2*C3-C1*S3 DER(2,5)=-C1*S2*C3 DER(3,5)=-C1*C2*S3-S1*C3 DER(1,6)=0. DER(2,6)=C2*C3 DER(3,6)=-S2*S3 DER(1,7)=C1*S2 DER(2,7)=S1*C2 DER(3,7)=0. DER(1,8)=S1*S2 DER(2,8)=-C1*C2 DER(3,8)=0. DER(1,9)=0. DER(2,9)=-S2 DER(3,9)=0. RETURN C END SUBROUTINE REFEUL (XT,XM,R,THE,NP) DIMENSION XT(3,30),XM(3,30),R(3,4),THE(3),AA(6,6 1), VEC(6), XDASH(3), DXYZ(3,6), DER(3,9), RCALC(9), SHIFT(6), EUL( 23), DEUL(3) WRITE (6,118) CONST=57.288 SIGOLD=0.0 DO 115 ICY=1,10 CALL CALC (THE,RCALC,DER,1) K=0 DO 101 I=1,3 DO 101 J=1,3 K=K+1 R(I,J)=RCALC(K) 101 CONTINUE DO 102 I=1,6 VEC(I)=0.0 DO 102 J=1,6 AA(I,J)=0.0 102 CONTINUE SIGMA=0.0 NCOUNT=0 P=1.0 DO 110 N=1,NP DO 103 I=1,3 DO 103 J=1,6 DXYZ(I,J)=0.0 IF (I.EQ.(J-3)) DXYZ(I,J)=1.0 103 CONTINUE DO 106 I=1,3 LB=(I-1)*3 XDASH(I)=R(I,4) DO 105 J=1,3 XDASH(I)=XDASH(I)+XM(J,N)*R(I,J) DO 104 K=1,3 L=LB+K DXYZ(I,J)=DXYZ(I,J)+DER(J,L)*XM(K,N) 104 CONTINUE 105 CONTINUE 106 CONTINUE DO 109 L=1,3 DIFF=XT(L,N)-XDASH(L) DO 108 I=1,6 VEC(I)=VEC(I)+DIFF*DXYZ(L,I)*P DO 107 J=1,6 AA(I,J)=AA(I,J)+DXYZ(L,I)*DXYZ(L,J)*P 107 CONTINUE 108 CONTINUE SIGMA=SIGMA+DIFF*DIFF 109 CONTINUE NCOUNT=NCOUNT+1 110 CONTINUE CALL SIMINV (AA,6) DO 112 I=1,6 SHIFT(I)=0.0 DO 111 J=1,6 SHIFT(I)=SHIFT(I)+AA(I,J)*VEC(J) 111 CONTINUE 112 CONTINUE DO 113 I=1,3 EUL(I)=THE(I)*CONST DEUL(I)=SHIFT(I)*CONST THE(I)=THE(I)+SHIFT(I) R(I,4)=R(I,4)+SHIFT(I+3) 113 CONTINUE RMS=SQRT(SIGMA/NCOUNT) WRITE (6,119) ICY,EUL,DEUL,(SHIFT(I),I=4,6),RMS IF (SIGOLD.EQ.0) GO TO 114 IF (((SIGOLD-SIGMA)/SIGOLD).LT.0.001) GO TO 116 114 SIGOLD=SIGMA 115 CONTINUE 116 CALL CALC (THE,RCALC,DER,0) K=0 DO 117 I=1,3 DO 117 J=1,3 K=K+1 R(I,J)=RCALC(K) 117 CONTINUE NCY=ICY RETURN C 118 FORMAT ( 48H REFINE EULERIAN ANGLES TO GIVE MIN FOR SUM DSQ ) 119 FORMAT ( 7H0CYCLE ,I2/, 17H EULERIAN ANGLES ,3F8.2/, 8H SHIFTS , 13F8.2/, 22H TRANSLATIONAL SHIFTS ,3F8.2/, 5H RMS ,F8.2) C END SUBROUTINE ROTANG (R,THE) DIMENSION R(3,4), THE(3), DER(3,9), ROBS(9), RCALC(9), AA(3,3), X( 13), DEL(3), EUL(3), DEUL(3) C C FIND APPROXIMATE VALUES OF THETA1,2,3. THETA2 LIES BETWEEN 0 AND 1 C WRITE (6,118) CONST=57.288 C2=R(3,3) IF (C2.GT.1.0) C2=1.0 IF (C2.LT.-1.0) C2=-1.0 THE(2)=ACOS(C2) IF ((1.0-ABS(C2)).LT.0.001) GO TO 101 S2=SIN(THE(2)) C2=COS(THE(2)) S1=R(3,1)/S2 C1=-R(3,2)/S2 S3=R(1,3)/S2 C3=R(2,3)/S2 THE(1)=ATAN2(S1,C1) THE(3)=ATAN2(S3,C3) GO TO 104 101 C11=R(1,1) IF (C11.GT.1.0) C11=1.0 IF (C11.LT.-1.0) C11=-1.0 C12=R(1,2) IF (C12.GT.1.0) C12=1.0 IF (C12.LT.-1.0) C12=-1.0 IF (C2.LT.0) GO TO 102 SUM13=ACOS(C11) DIF13=ACOS(C12) GO TO 103 102 SUM13=ACOS(-C12) DIF13=ACOS(C11) 103 THE(1)=(SUM13+DIF13)/2.0 THE(3)=(SUM13-DIF13)/2.0 104 K=0 DO 106 I=1,3 DO 105 J=1,3 K=K+1 ROBS(K)=R(I,J) 105 CONTINUE 106 CONTINUE DO 114 ICY=1,15 DO 108 I=1,3 X(I)=0.0 DO 107 J=1,3 AA(I,J)=0.0 107 CONTINUE 108 CONTINUE CALL CALC (THE,RCALC,DER,1) DO 111 K=1,9 DO 110 I=1,3 X(I)=X(I)+(ROBS(K)-RCALC(K))*DER(I,K) DO 109 J=1,3 AA(I,J)=AA(I,J)+DER(I,K)*DER(J,K) 109 CONTINUE 110 CONTINUE 111 CONTINUE CALL SIMINV (AA,3) SUM=0.0 DO 113 I=1,3 DEL(I)=0.0 DO 112 J=1,3 DEL(I)=DEL(I)+AA(I,J)*X(J) 112 CONTINUE EUL(I)=THE(I)*CONST THE(I)=THE(I)+DEL(I) DEUL(I)=DEL(I)*CONST SUM=SUM+ABS(DEL(I)) 113 CONTINUE WRITE (6,119) ICY,EUL,DEUL,SUM IF (SUM.LT.0.003) GO TO 115 114 CONTINUE ICY=15 115 CALL CALC (THE,RCALC,DER,0) SUM1=0.0 SUM2=0.0 K=0 DO 117 I=1,3 DO 116 J=1,3 K=K+1 SUM1=SUM1+ABS(ROBS(K)-RCALC(K)) SUM2=SUM2+ABS(ROBS(K)) R(I,J)=RCALC(K) 116 CONTINUE 117 CONTINUE SUM1=SUM1/SUM2 WRITE (6,120) ICY,SUM1 RETURN C 118 FORMAT ( 59H FIND BEST EULERIAN ANGLES FROM LINEARLY DETERMINED MA 1TRIX ) 119 FORMAT ( 7H CYCLE ,I2, 17H EULERIAN ANGLES ,3F8.2, 7HSHIFTS ,3F8 1.2, 5H SUM ,F8.2) 120 FORMAT (1X,I2, 25H CYCLES, FINAL RESIDUAL ,F8.4) C END SUBROUTINE SIMINV (B,N) C C INVERTS REAL SYMMETRIC MATRIX UPTO ORDER 6. C LOWER TRIANGLE IN ARRAY *A* C DIMENSION B(N,N) DIMENSION A(21), P(6), Q(6), R(6) LOGICAL R M=1 DO 101 I=1,N DO 101 J=1,I A(M)=B(I,J) M=M+1 101 CONTINUE DO 102 I=1,N 102 R(I)=.TRUE. DO 112 I=1,N BIGAJJ=0.0 KK=0 DO 104 J=1,N KK=KK+J IF (R(J).AND.(ABS(A(KK)).GT.BIGAJJ)) GO TO 103 GO TO 104 103 BIGAJJ=ABS(A(KK)) K=J JJ=KK 104 CONTINUE IF (BIGAJJ.NE.0.0) GO TO 105 PRINT 114 Q(K)=0.0 GO TO 106 105 Q(K)=1.0/A(JJ) 106 R(K)=.FALSE. P(K)=1.0 A(JJ)=0.0 IF (K.EQ.1) GO TO 108 M=K-1 JK=JJ-K DO 107 J=1,M JK=JK+1 P(J)=A(JK) Q(J)=A(JK)*Q(K) IF (R(J)) Q(J)=-Q(J) A(JK)=0.0 107 CONTINUE IF (K.EQ.N) GO TO 110 108 M=K+1 JK=JJ DO 109 J=M,N JK=JK+J-1 P(J)=-A(JK) IF (R(J)) P(J)=-P(J) Q(J)=-A(JK)*Q(K) A(JK)=0.0 109 CONTINUE 110 JK=0 DO 111 J=1,N DO 111 K=1,J JK=JK+1 111 A(JK)=A(JK)+P(K)*Q(J) 112 CONTINUE M=1 DO 113 I=1,N DO 113 J=1,I B(I,J)=A(M) B(J,I)=A(M) M=M+1 113 CONTINUE RETURN C 114 FORMAT (14H SYMINV FAILED) C END