PROGRAM REDUCE C C BUERGER CELL REDUCTION C C L. V. AZAROFF AND M. J. BUERGER, THE POWDER METHOD IN X-RAY C CRYSTALLOGRAPHY, MC GRAW-HILL, NEW YORK, 1958, PP. 125-159. C DIMENSION A(3),ALPHA(3),U(3),G(3,3),TT(0:2,-2:+2,-2:+2) DIMENSION T(3,3),V(3),B(3),BETA(3) DIMENSION DATA(75),INDEX(75) PI=3.141593 C C START WITH AN ARBITRARY PRIMITIVE CELL. C --------- C IO1=10 IO5=5 IO6=6 OPEN(UNIT=IO1,FILE='cell.dat',STATUS='UNKNOWN') WRITE (IO6,999) 999 FORMAT (/1X,'CELL PARAMETERS FROM "cell.dat" FILE? Y OR N.') READ (IO5,899) ANS 899 FORMAT (A) 100 IF (ANS.NE.'Y') THEN WRITE (IO6,896) READ (IO5,*) A,ALPHA REWIND IO1 WRITE (IO1,895) A,ALPHA ELSE READ (IO1,*) A,ALPHA END IF 896 FORMAT (/1X,'ENTER A, B, C, ALPHA, BETA, GAMMA.') 895 FORMAT (3F10.5,3F10.4) CA=COS(ALPHA(1)*PI/180) CB=COS(ALPHA(2)*PI/180) CG=COS(ALPHA(3)*PI/180) VOL=A(1)*A(2)*A(3)*(1-CA**2-CB**2-CG**2+2*CA*CB*CG) WRITE (IO6,898) A,ALPHA,VOL 898 FORMAT (/1X,' A B C ALPHA BETA ', &' GAMMA V'/1X,3F10.4,3F10.3,F10.2) WRITE (IO6,897) 897 FORMAT (1X,'CORRECT? Y OR N.') READ (IO5,899) ANS IF (ANS.NE.'Y') GO TO 100 C C STORE SCALAR PRODUCTS OF LATTICE TRANSLATION VECTORS, I.E., THE C METRIC TENSOR G(I,J). C G(1,1)=A(1)**2 G(1,2)=A(1)*A(2)*COS(ALPHA(3)*PI/180) G(1,3)=A(1)*A(3)*COS(ALPHA(2)*PI/180) G(2,1)=G(1,2) G(2,2)=A(2)**2 G(2,3)=A(2)*A(3)*COS(ALPHA(1)*PI/180) G(3,1)=G(1,3) G(3,2)=G(2,3) G(3,3)=A(3)**2 C C STORE IN MATRIX FORM TABLES OF T*T VALUES, WHERE C T = U*A + V*B + W*C, FOR U = 0, 1, 2 AND V, W = -2, -1, 0, 1, 2. C DO 20 IU= 0,+2 DO 20 IV=-2,+2 DO 20 IW=-2,+2 20 TT(IU,IV,IW)=0 DO 21 IU= 0,+2 U(1)=IU DO 21 IV=-2,+2 U(2)=IV DO 21 IW=-2,+2 U(3)=IW C C FUNCTION VTMV CALCULATES THE BILINEAR FORM, C SCALAR = (VECTOR TRANSPOSE)*(MATRIX)*(VECTOR). C 21 TT(IU,IV,IW)=VTMV(U,G,U) C C PRINT THE T*T TABLES. C ILP=60 OPEN (UNIT=ILP,FILE='reduce.lp',STATUS='NEW') WRITE (ILP,997) (A(I),ALPHA(I),I=1,3),VOL 997 FORMAT (1H1,'CELL REDUCTION'//1X,'A = ',F10.4,' ALPHA = ', &F10.3/' B = ',F10.4,' BETA = ',F10.3/' C = ',F10.4,' GAMM', &'A = ',F10.3//' V = ',F10.2) WRITE (ILP,996) 996 FORMAT (/1X,'MATRIX TABLES OF VALUES OF T*T, WHERE T = U*A + V*B + & W*C.'/1X,'TO FIND THE THREE SHORTEST NONCOPLANAR LATTICE TRANSLAT &ION VECTORS,'/' LOOK FOR THE U, V, W FOR THE THREE SMALLEST T*T VA &LUES.') DO 22 IU=0,2 WRITE (ILP,995) IU 995 FORMAT (/1X,'U = ',I2,8X,'V = -2',8X,'-1',9X,'0',9X,'1',9X,'2') IW=-2 WRITE (ILP,994) IW,(TT(IU,IV,IW),IV=-2,+2) 994 FORMAT (/1X,'W = ',I2,4X,5F10.3) DO 22 IW=-1,+2 WRITE (ILP,993) IW,(TT(IU,IV,IW),IV=-2,+2) 993 FORMAT (5X,I2,4X,5F10.3) 22 CONTINUE C C PRINT SORTED LIST OF THE T*T MATRIX ELEMENTS. C WRITE (IO6,'(/1X,''PRINT SORTED LIST OF THE T*T MATRIX ELEMENTS'', &''? Y OR N.'')') READ (IO5,'(A1)') ANS IF (ANS.EQ.'Y') THEN WRITE (ILP,970) N=0 DO IU= 0,+2 DO IV=-2,+2 DO IW=-2,+2 N=N+1 DATA(N)=TT(IU,IV,IW) END DO END DO END DO CALL SORT (N,DATA,INDEX) DO I=1,N DO IU= 0,+2 DO IV=-2,+2 DO IW=-2,+2 IF (DATA(INDEX(I)).EQ.TT(IU,IV,IW)) WRITE (ILP,971) & IU,IV,IW,TT(IU,IV,IW) END DO END DO END DO END DO END IF 970 FORMAT (/1X,' SORTED LIST OF T*T MATRIX ELEMENTS:'//1X,' IU IV I', &'W T*T(IU,IV,IW)'/1X,' -- -- -- -------------') 971 FORMAT (1X,3I3,F10.3) C C TRANSFORM TO NEW UNIT CELL. C WRITE (IO6,980) 980 FORMAT (/1X,'ENTER THE TRANSFORMATION MATRIX? Y OR N.') READ (IO5,981) ANS 981 FORMAT (A1) IF (ANS.NE.'Y') STOP WRITE (IO6,992) 992 FORMAT (/1X,'NEW A = T11*A + T12*B + T13*C'/' NEW B = T21*A + T22* &B + T23*C'/' NEW C = T31*A + T32*B + T33*C'//1X,'ENTER T11, T12, T &13.') READ (IO5,*) T(1,1),T(1,2),T(1,3) WRITE (IO6,991) 991 FORMAT (' ENTER T21, T22, T23.') READ (IO5,*) T(2,1),T(2,2),T(2,3) WRITE (IO6,990) 990 FORMAT (' ENTER T31, T32, T33.') READ (IO5,*) T(3,1),T(3,2),T(3,3) WRITE (ILP,989) ((T(I,J),J=1,3),I=1,3) 989 FORMAT (/1X,'TRANSFORMATION FROM GIVEN TO TRANSFORMED LATTICE TRAN &SLATION VECTORS'//1X,'T11 T12 T13 ',3F5.1/1X,'T21 T22 T23 ', &3F5.1/' T31 T32 T33 ',3F5.1) DO 31 I=1,3 DO 32 J=1,3 32 U(J)=T(I,J) 31 B(I)=SQRT(VTMV(U,G,U)) DO 41 I=1,3 J=I+1 K=I+2 IF (J.GT.3) J=J-3 IF (K.GT.3) K=K-3 DO 42 L=1,3 U(L)=T(J,L) 42 V(L)=T(K,L) 41 BETA(I)=ACOS(VTMV(U,G,V)/SQRT(VTMV(U,G,U)*VTMV(V,G,V)))*180/PI CA=COS(BETA(1)*PI/180) CB=COS(BETA(2)*PI/180) CG=COS(BETA(3)*PI/180) VOL=B(1)*B(2)*B(3)*(1-CA**2-CB**2-CG**2+2*CA*CB*CG) WRITE (ILP,988) (B(I),BETA(I),I=1,3),VOL 988 FORMAT (/1X,'TRANSFORMED CELL PARAMETERS'//1X,'A = ',F10.4,' AL &PHA = ',F10.3/1x,'B = ',F10.4,' BETA = ',F10.3/1x,'C = ', &F10.4,' GAMMA = 'F10.3//1X,'V = ',F10.2) WRITE (IO1,895) B,BETA CLOSE (UNIT=IO1) C C EVALUATE AND PRINT NIGGLI MATRIX OF SCALAR PRODUCTS. C S11=B(1)**2 S22=B(2)**2 S33=B(3)**2 S23=B(2)*B(3)*COS(BETA(1)*PI/180) S31=B(3)*B(1)*COS(BETA(2)*PI/180) S12=B(1)*B(2)*COS(BETA(3)*PI/180) WRITE (ILP,987) S11,S22,S33,S23,S31,S12 987 FORMAT (/1X,'NIGGLI MATRIX FOR THE TRANSFORMED CELL'//1X,'S11 S2', &'2 S33 ',3F10.4/1X,'S23 S31 S12 ',3F10.4) STOP END C----------------------------------------------------------------------- FUNCTION VTMV(U,G,V) C C BILINEAR FORM = (VECTOR TRANSPOSE)*(MATRIX)*(VECTOR) = SCALAR C BILINEAR FORM = (ROW MATRIX)*(SQUARE MATRIX)*(COLUMN MATRIX) C DIMENSION U(3),G(3,3),V(3) VTMV=0 DO 1 I=1,3 DO 1 J=1,3 1 VTMV=VTMV+U(I)*G(I,J)*V(J) RETURN END C------------------------------------------------------------------------ SUBROUTINE SORT (N,DATA,INDEX) C C INDEXES THE ARRAY DATA(N) AND RETURNS THE ARRAY INDEX(N) SORTED SUCH C THAT THE VALUES DATA(INDEX(I)) ARE IN ASCENDING ORDER FOR I = 1, 2, C ..., N. THE INPUT VARIABLES N AND DATA(N) ARE NOT CHANGED. C C EMPLOYS THE "HEAPSORT" ALGORITHM. C C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A. C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986). NUMERICAL RECIPIES: THE C --------- -------- --- C ART OF SCIENTIFIC COMPUTING, PP. 229-233. CAMBRIDGE, ENGLAND: C --- -- ---------- --------- C CAMBRIDGE UNIVERSITY PRESS. C REAL DATA,T DIMENSION DATA(N),INDEX(N) DO I=1,N INDEX(I)=I END DO IF (N.EQ.1) RETURN I=N/2+1 J=N 1 CONTINUE IF (I.GT.1) THEN I=I-1 INEXT=INDEX(I) T=DATA(INEXT) ELSE INEXT=INDEX(J) T=DATA(INEXT) INDEX(J)=INDEX(1) J=J-1 IF (J.EQ.1) THEN INDEX(1)=INEXT RETURN END IF END IF K=I L=I+I DO WHILE (L.LE.J) IF (L.LT.J) THEN IF (DATA(INDEX(L)).LT.DATA(INDEX(L+1))) L=L+1 END IF IF (T.LT.DATA(INDEX(L))) THEN INDEX(K)=INDEX(L) K=L L=L+L ELSE L=J+1 END IF END DO INDEX(K)=INEXT GO TO 1 END C-----------------------------------------------------------------------