C----------------------------------------------------------------------- C This program finds the conventional representation of a lattice C input as cell parameters and a lattice type, assuming that metric C relations in the lattice correspond to lattice symmetry. C Pseudo-symmetry in the primitive lattice is also detected. C See: Le Page, Y. (1982). J. Appl. Cryst., 15,255-259. C Sept. 1986 Fortran 77 + three-fold axes YLP. C----------------------------------------------------------------------- SUBROUTINE CREDUC (KI) COMMON /GEOM/ AA(3,3),AINV(3,3),TRANS(3,3),RH(3,20),HH(3,20), $ AANG(20),PH(3,20),PMESH(3,2,20),NERPAX(20),N2,N3, $ EXPER COMMON /IOUASS/ IOUNIT(10) CHARACTER COUT*132 COMMON /IOUASC/ COUT(20) COMMON /IODEVS/ ITP,ITR,LPT,LPTX,LNCNT,PGCNT,ICD,IRE,IBYLEN, $ IPR,NPR,IIP COMMON /INFREE/ IFREE(20),RFREE(20),ICFLAG DIMENSION P(3),H(3),IP(3),IR(3),HX(3,37),DHX(3,37),CELL(3), $ CELANG(3),SHORT(4,4),IPAD(20),VPROD(3),DIRECT(3), $ RECIP(3) CHARACTER KI*2 C---------------------------------------------------------------------- C The 37 acceptable index combinations of 0, +/- 1 or 2 C---------------------------------------------------------------------- DATA ITOT/37/ DATA DHX/ 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, $ 0, 1, 1, 1,-1, 0, 1, 0,-1, 0, 1,-1, 1, 1, 1, $ 1, 1,-1, 1,-1, 1, -1, 1, 1, 2, 1, 0, 2, 0, 1, $ 2,-1, 0, 2, 0,-1, 0, 2, 1, 1, 2, 0, 0, 2,-1, $ -1, 2, 0, 1, 0, 2, 0, 1, 2, -1, 0, 2, 0,-1, 2, $ 2, 1, 1, 2, 1,-1, 2,-1, 1, -2, 1, 1, 1, 2, 1, $ 1, 2,-1, 1,-2, 1, -1, 2, 1, 1, 1, 2, 1, 1,-2, $ 1,-1, 2, -1, 1, 2/ INP = -1 IOUT = ITP WRITE (COUT,18000) CALL FREEFM (ITR) EXPER = RFREE(1) IF (EXPER .EQ. 0.0) EXPER = 0.01 C---------------------------------------------------------------------- C Get the input cell and bring back the Buerger cell parameters C and the input -> Buerger cell parameters C---------------------------------------------------------------------- 100 CALL CINPUT (IOUT,CELL,CELANG,TRANS) WRITE (COUT,13000) CALL GWRITE (IOUT,' ') WRITE (COUT,14000) CALL GWRITE (IOUT,' ') C----------------------------------------------------------------------- C Describe the Buerger direct and reciprocal cells by their cartesian C coordinates AA and AINV. CRAP is a dummy floating argument C----------------------------------------------------------------------- CALL MATRIX (CELL,CELANG,AA,CRAP,'ORTHOG') CALL MATRIX (AA,AINV,CRAP,CRAP,'INVERT') C---------------------------------------------------------------------- C Default angular tolerance: 3 degrees C---------------------------------------------------------------------- ANGMAX = TAN(3.0/57.2958)**2 C----------------------------------------------------------------------- C Find the twofold axes: C Generate all unique combinations of 0, 1 and 2 C Get the direction cosines of the possible rows C----------------------------------------------------------------------- DO 110 IT = 1,ITOT CALL MATRIX (AA,DHX(1,IT),HX(1,IT),CRAP,'MATVEC') 110 CONTINUE C----------------------------------------------------------------------- C Get direction cosines for the normal to the possible planes in turn C----------------------------------------------------------------------- N2 = 0 DO 140 IT = 1,ITOT CALL MATRIX (DHX(1,IT),AINV,P,CRAP,'VECMAT') C----------------------------------------------------------------------- C Select the rows in turn C----------------------------------------------------------------------- DO 130 L = 1,ITOT C----------------------------------------------------------------------- C Calculate the multiplicity of the cell defined by the mesh on the C plane and the translation along the row C----------------------------------------------------------------------- MULT = ABS(DHX(1,L)*DHX(1,IT) + DHX(2,L)*DHX(2,IT) + $ DHX(3,L)*DHX(3,IT)) + 0.1 IF (MULT .EQ. 1 .OR. MULT .EQ. 2) THEN C----------------------------------------------------------------------- C Calculate the angle between the row and the normal to the plane C----------------------------------------------------------------------- ANG = ((P(1)*HX(2,L) - P(2)*HX(1,L))**2 + $ (P(2)*HX(3,L) - P(3)*HX(2,L))**2 + $ (P(3)*HX(1,L) - P(1)*HX(3,L))**2) IF (ANG .LE. ANGMAX) THEN N2 = N2 + 1 DO 120 NX = 1,3 PH(NX,N2) = DHX(NX,IT) HH(NX,N2) = HX (NX,L) RH(NX,N2) = DHX(NX,L) 120 CONTINUE AANG(N2) = ANG NERPAX(N2) = 2 ENDIF ENDIF 130 CONTINUE 140 CONTINUE N3 = N2 C----------------------------------------------------------------------- C Order the rows on the angle with the normal to the plane C----------------------------------------------------------------------- IF (N2 .LT. 2) GO TO 250 DO 170 I = 1,N2 - 1 ANMAX = AANG(I) MAX = I DO 160 J = I + 1,N2 IF (AANG(J) .LT. ANMAX) THEN ANMAX = AANG(J) MAX = J ENDIF 160 CONTINUE CALL MATRIX (RH(1,I),RH(1,MAX),CRAP,CRAP,'INTRCH') CALL MATRIX (PH(1,I),PH(1,MAX),CRAP,CRAP,'INTRCH') CALL MATRIX (HH(1,I),HH(1,MAX),CRAP,CRAP,'INTRCH') AANG(MAX) = AANG(I) AANG(I) = ANMAX 170 CONTINUE C----------------------------------------------------------------------- C Find the kind of axis: find a family of coplanar twofold axes C----------------------------------------------------------------------- IF (N2 .LT. 3) GO TO 250 DO 220 I = 1,N2 - 1 IPAD(1) = I DO 210 J = I + 1,N2 IPAD(2) = J NUMAX = 2 DO 180 K = 1,N2 IF (K .NE. I .AND. K .NE. J) THEN CALL MATRIX (RH(1,I),RH(1,J),RH(1,K),DET,'DETERM') IF (ABS(DET) .LE. 0.01) THEN IF (K .LT. J) GO TO 210 NUMAX = NUMAX + 1 IPAD (NUMAX) = K ENDIF ENDIF 180 CONTINUE C----------------------------------------------------------------------- C Now find a twofold axis perpendicular to this plane C----------------------------------------------------------------------- DO 190 K = 1,N2 CALL MATRIX (PH(1,K),RH(1,I),SCAL,CRAP,'SCALPR') IF (ABS(SCAL) .LE. 0.01) THEN CALL MATRIX (PH(1,K),RH(1,J),SCAL,CRAP,'SCALPR') IF (ABS(SCAL) .LE. 0.01) THEN C----------------------------------------------------------------------- C Found one: its maximum order is NUMAX, the number of perpend. axes C----------------------------------------------------------------------- NERPAX(K) = NUMAX IF (NUMAX .LE. 2) NERPAX(K) = 2 GO TO 210 ENDIF ENDIF 190 CONTINUE C----------------------------------------------------------------------- C Three coplanar axes were found, but no perpendicular one: C this is likely to be a threefold axis. C----------------------------------------------------------------------- IF (NUMAX .GT. 2) THEN CALL MATRIX (HH(1,I),HH(1,J),VPROD ,CRAP,'VECPRD') CALL MATRIX (VPROD ,AA ,RECIP ,CRAP,'VECMAT') CALL MATRIX (RECIP ,RECIP ,CRAP ,CRAP,'COPRIM') CALL MATRIX (RECIP ,AINV ,P ,CRAP,'VECMAT') CALL MATRIX (AINV ,VPROD ,DIRECT,CRAP,'MATVEC') CALL MATRIX (DIRECT ,DIRECT ,CRAP ,CRAP,'COPRIM') CALL MATRIX (AA ,DIRECT ,H ,CRAP,'MATVEC') CALL MATRIX (DIRECT ,RECIP ,SCAL ,CRAP,'SCALPR') MULT = ABS(SCAL) + 0.1 ANG = ((P(1)*H(2) - P(2)*H(1))**2 + $ (P(2)*H(3) - P(3)*H(2))**2 + $ (P(3)*H(1) - P(1)*H(3))**2)/(MULT*MULT) IF (ANG .LE. ANGMAX) THEN C----------------------------------------------------------------------- C All seems to be ok, save the results C----------------------------------------------------------------------- N3 = N3 + 1 DO 200 NX = 1,3 PH(NX,N3) = RECIP(NX) RH(NX,N3) = DIRECT(NX) HH(NX,N3) = H(NX) 200 CONTINUE AANG(N3) = ANG NERPAX(N3) = NUMAX IF (NUMAX .EQ. 0) NERPAX(N3) = 2 ENDIF ENDIF 210 CONTINUE 220 CONTINUE C----------------------------------------------------------------------- C Order the threefold axes on the angle with the plane C----------------------------------------------------------------------- IF (N3 - N2 .GE. 2) THEN DO 240 I = N3,N3 - 1,-1 ANMAX = AANG(I) MAX = I DO 230 J = I + 1,N3 IF (AANG(J) .LT. ANMAX) THEN ANMAX = AANG(J) MAX = J ENDIF 230 CONTINUE CALL MATRIX (RH(1,I),RH(1,MAX),CRAP,CRAP,'INTRCH') CALL MATRIX (PH(1,I),PH(1,MAX),CRAP,CRAP,'INTRCH') CALL MATRIX (HH(1,I),HH(1,MAX),CRAP,CRAP,'INTRCH') SAVE = NERPAX(I) NERPAX(I) = NERPAX(MAX) NERPAX(MAX) = NERPAX(I) AANG(MAX) = AANG(I) AANG(I) = ANMAX 240 CONTINUE ENDIF C----------------------------------------------------------------------- C Get 2 primitive translations for the perpendicular plane C----------------------------------------------------------------------- 250 DO 380 IT = 1,N3 NMESH = 1 DO 260 I = 1,ITOT CALL MATRIX (DHX(1,I),PH(1,IT),SCAL,CRAP,'SCALPR') IF (ABS(SCAL) .LE. 0.01) THEN NMESH2 = I IF (NMESH .EQ. 1) NMESH1 = I IF (NMESH .EQ. 2) GO TO 270 NMESH = NMESH + 1 ENDIF 260 CONTINUE 270 DO 280 I = 1,3 SHORT(I,1) = DHX(I,NMESH1) SHORT(I,2) = DHX(I,NMESH2) 280 CONTINUE C----------------------------------------------------------------------- C Get the 2 shortest translations in the plane: generate mesh diagonals C----------------------------------------------------------------------- 290 DO 300 I = 1,3 SHORT(I,3) = SHORT(I,1) + SHORT(I,2) SHORT(I,4) = SHORT(I,1) - SHORT(I,2) 300 CONTINUE DO 310 I = 1,4 CALL MATRIX (AA,SHORT(1,I),SHORT(4,I),CRAP,'LENGTH') 310 CONTINUE C----------------------------------------------------------------------- C Rank their lengths C----------------------------------------------------------------------- ISWTCH = 0 DO 340 I = 1,2 DO 330 J = 2,4 IF (SHORT(4,J) .LT. SHORT(4,I)) THEN DO 320 K = 1,4 SAVE = SHORT(K,I) SHORT(K,I) = SHORT(K,J) SHORT(K,J) = SAVE 320 CONTINUE ISWTCH = 1 ENDIF 330 CONTINUE 340 CONTINUE C----------------------------------------------------------------------- C Finished when no more interchanges C----------------------------------------------------------------------- IF (ISWTCH .EQ. 1) GO TO 290 C----------------------------------------------------------------------- C Make sure the angle is not acute C----------------------------------------------------------------------- CALL MATRIX (AA,SHORT(1,1),SHORT(1,3),CRAP,'MATVEC') CALL MATRIX (AA,SHORT(1,2),SHORT(1,4),CRAP,'MATVEC') CALL MATRIX (SHORT(1,3),SHORT(1,4),SCAL,CRAP,'SCALPR') IF (SCAL .GE. 0.0) THEN DO 350 IAX = 1,3 SHORT(IAX,2) = -SHORT(IAX,2) 350 CONTINUE ENDIF C----------------------------------------------------------------------- C Make sure the reference system is right-handed C----------------------------------------------------------------------- IS = 1 CALL MATRIX (RH(1,IT),PH(1,IT),SCAL,CRAP,'SCALPR') IF (SCAL .LT. 0.) IS = -1 IS1 = 1 CALL MATRIX (SHORT(1,1),SHORT(1,2),RH(1,IT),DET,'DETERM') IF (DET .LT. 0.) IS1 = -1 C----------------------------------------------------------------------- C This is a potential symmetry axis, we print and save the values C----------------------------------------------------------------------- DO 370 NX = 1,3 RH(NX,IT) = IS1*RH(NX,IT) HH(NX,IT) = IS1*HH(NX,IT) PH(NX,IT) = IS*IS1*PH(NX,IT) IP(NX) = PH(NX,IT) IR(NX) = RH(NX,IT) DO 370 NY = 1,2 PMESH (NX,NY,IT) = IS1 * SHORT (NX,NY) 370 CONTINUE AANG(IT) = ATAN(SQRT(AANG(IT)))*180.0/3.1415927 MULT = IR(1)*IP(1) + IR(2)*IP(2) + IR(3)*IP(3) WRITE (COUT,15000) IR,IP,MULT,AANG(IT),NERPAX(IT) CALL GWRITE (IOUT,' ') 380 CONTINUE C----------------------------------------------------------------------- C Fill the next slot C----------------------------------------------------------------------- DO 390 I = 1,3 RH(I,N3 + 1) = 0.0 PH(I,N3 + 1) = 0.0 PMESH(I,1,N3 + 1) = 0.0 PMESH(I,2,N3 + 1) = 0.0 390 CONTINUE PMESH(1,1,N3 + 1) = 1.0 PMESH(2,2,N3 + 1) = 1.0 RH(3,N3 + 1) = 1.0 PH(3,N3 + 1) = 1.0 C----------------------------------------------------------------------- C Find the crystal system C----------------------------------------------------------------------- WRITE (COUT,16000) CALL GWRITE (IOUT,' ') NPSUDO = N2 CALL FNDSYS (IOUT,HH,NPSUDO) IF (INP .GT. 0 .OR. IOUT .NE. ITP) THEN WRITE (COUT,17000) CALL GWRITE (IOUT,' ') ENDIF KI = ' ' RETURN 9000 FORMAT (/10X,'CREDUC -- The NRCVAX Cell Reduction Routine'/'%') 10000 FORMAT (' Input from the terminal or a file (T) ? ',$) 11000 FORMAT (' Output to terminal or lineprinter-file (T) ? ',$) 13000 FORMAT (/15X,'Possible 2-fold Axes:'/ $ 14X,'Rows',20X,'Products',9X,'Kind') 14000 FORMAT (7X,'Direct',6X,'Reciprocal',7X,'Dot',4X,'Vector',4X, $ 'of Axis') 15000 FORMAT (2X,3I4,2X,3I4,I10,F10.3,7X,I3) 16000 FORMAT (/) 17000 FORMAT (//) 18000 FORMAT (' Type the Allowable Tolerance on True Cell Angles', $ ' (0.01deg) ',$) END