C
C     FILE    SPELIB                                Last Update: 1986.12.11
C             ------                                            I.D.Brown
C
C     THIS FILE CONTAINS:
C
C  S SPOST(IER,ND)            DERIVES ONE CORDINATE TRIPLET FOR EACH SPECIAL
C                             POSITION
C  S REDEFI(SPS,TPS)          REDEFINES INDEPENDENT VARIABLES IN SEITZ MATRIS
C                             (SPS/TPS) TO INCLUDE TRANSLATIONSc
C  S NEWAXS(SPS,TPS)          REASSIGNES AXES OF INDEPENDENT VARIABLES IN
C                             SEITZ MATRIX
C  S TEHECU(SP,SPS)           REASSIGNES DEPENDENT AXES IN SEITZ MATRIX FOR
C                             HIGHER SYMMETRIES
C  S SHCTOR(TP,SP,TU,NLT)     REDEFINES TRANSLATIONS IN SEITZ MATRIX WITH
C                             RESPECT TO LATTICE TRANSLATIONS
C  S SPMULT                   CALCULATES THE MULTIPLICATION TABLE FOR SPECIAL
C                             POSITIONS
C  S SPESOR                   DOES STANDARD ORDERING OF SPECIAL POSITIONS
C  S SPSYM                    CREATES POINT GROUP SYMBOL FOR ALL SPECIAL
C                             POSITIONS
C  F ITESTR(TA,TB,TU,NLT,SP)  CHECKS ON EQUIVALENCE OF TWO TRANSLATIONS
C  S RECENT(X,N)              BRINGS N-VECTOR ELEMENTS INTO RANGE (.0,.99)
C  S LTRANS(NLT,T,TU)         ADDS CENTRING TRANSLATIONS TO ARRAY =TU=
C  S XALTNO(XX,NPOS,NA1,NA2,IC,NT,ZZ) ZZ IS THE SAME AS XX TRANFORMED TO
C                                     THE ALTERMATT POSITION
C  S NALTNO(XX,NOS,NA1,NA2,IC,NT)   DERIVES THE SPECIAL POSITION NUMBER
C                                   FOR AN ATOMIC POSITION
C  F IALTST(SN,TM,X)          "SLAVE" TO NALTNO
C
C     EXTERNALS CALLED IN THIS FILE:
C
C  MATHLB: ICOMMA, ICOMVE, INTEGER, MMINM, SECOND, VMINV, VPLUSV
C  SYMLIB: SYMOP, SYMULT
C
C
C****** STANDARD FORTRAN 77 CODE ONLY USED IN THIS FILE
C
C
C
C
      SUBROUTINE SPOST(IER,NBUG)
C     --------------------------
C
C   WRITTEN BY D. ALTERMATT
C   MCMASTER UNIVERSITY, OCT/NOV 1984
C
C   THIS SUBROUTINE CALCULATES THE SPECIAL POSITION TABLE STORED IN
C   ISPTAB
C
C   CALLS LTRANS, VPLUSV, MMINM, SYMULT, ICOMMA, SHCTOR,
C         RECENT, ITESTR, INTEGER, ICOMVE, REDEFI, TEHECU,
C         NEWAXS
C
C***** CALLS TIME-FUNCTION =SECOND(I)=
C
C                    THEORY
C                    ------
C
C   IF THE VECTOR X REPRESENTS THE COORDINATES OF A SPECIAL POSITION
C   THEN
C          X = SX+T+L
C   WHERE S/T IS A SYMMETRY OPERATOR AND L IS ANY UNIT LATTICE TRANSLATION
C
C   HENCE
C          (I-S)X = T+L
C   OR          X = (I-S)**(-1)*(T+L)
C
C   EQUATIONS WITHOUT SOLUTIONS CORRESPOND TO GLIDES AND SCREWS
C   THAT DO NOT GIVE RISE TO SPECIAL POSITIONS.
C   IN MANY CASES EQUATIONS ARE INDETERMINATE BECAUSE A COLUMN 'I'
C   AND ROW 'I' ARE ALL ZERO OR MATRIX S IS SINGULAR. IN THESE CASES
C   THE VALUE OF X(I) IS NOT FIXED OR DEPENDENT ON X(J) BY SYMMETRY.
C
      CHARACTER*8 GLAT,GCENT,GCFL
C
C   SYSTEM COMMON BLOCKS
C
      COMMON/FILES/NIN,NOUT,NDEBUG
      COMMON/SYM/NSYM,MSYM,S(3,3,48),T(3,48),NTT,TA(3,4)
      COMMON/GSYM/GLAT,GCENT
      COMMON/SYMA/MULT(48,48),IT(3,48),NSP,MSP,ISPTAB(48,30),
     +   MULTSP(30),SP(3,3,30),TP(3,30)
C
      DIMENSION SD(3,3),TD(3),TU(3,14),LR(3),L2(3),LC(3),SS(3,3),
     + IPOINT(48),ROW(3),AROW(3),COL(3),ACOL(3),KZ(3),TT(3),TPP(3)
      DIMENSION ISAVE(48),SPS(3,3,48),TPS(3,48),R222(3),TS(3),TSS(3)
      DIMENSION SPP(3,3),KZZ(3)
C   GENERAL LATTICE TRANSLATIONS
      DATA TU/1.,1.,1., 1.,0.,0., 0.,1.,0., 0.,0.,1., 0.,1.,1.,
     + 1.,0.,1., 1.,1.,0., 0.,0.,0., 18*.0/
C
C   TO AVOID EXESS OUTPUT, THE SWITCH FOR DEBUG-OUTPUT =NDEBUG= IS
C   TEMPORARELLY REDEFINED AT THIS POINT
C
      NDEB=NDEBUG
      NDEBUG=NBUG
C
      MSP=30
      IER=0
      NLT=8
C   FLAGS TO FIND DEPENDANT AXIS
      ICRFL1=0
      ICRFL2=0
      ICRFL3=0
C   FLAGS TO FIND HEXAGONAL AXIS
      IHEXE=0
      IHEXZ=0
      IHEXS=0
      M3=0
C   FLAGS TO FIND 11Z-SYMMETRY (AND RELATED..)
      I11Z=0
      I11ZS=0
      I11ZC=0
C   FLAGS TO FIND 222-SYMMETRY
      R222(1)=0.0
      R222(2)=0.0
      R222(3)=0.0
      I222=0
      I222L=0
      I222S=0
      I2LIM=1
C   ZERO COUNTERS
      NSP=0
      NSP1=0
      NSP2=0
      NSP3=0
      NSP4=0
      NSP5=0
      NSP6=0
C   FLAG TO FIND 1/4,1/4,3/4
      I444=0
C   FLAG FOR CENTRE OF SYMMETRY (ANYWHERE)
      GCFL='A'
      IF(GCENT.NE.'C') THEN
        L=NSYM/2+1
        DO 5 I=1,L
          TRACE=S(1,1,I)+S(2,2,I)+S(3,3,I)
          OTR=S(1,2,I)+S(1,3,I)+S(2,3,I)+S(2,1,I)+S(3,1,I)+S(3,2,I)
          IF(INTEGER(TRACE,1).EQ.-30.AND.INTEGER(OTR,3).EQ.0) THEN
            GCFL='C'
            GOTO 10
            ENDIF
    5     CONTINUE
      ELSE
        GCFL='C'
        ENDIF
C
C   FIRST BRING ALL TRANSLATIONS INTO RANGE 0 TO .99 AND RESET
C   SYM-OP POINTER (.NE.1 IF SYMMETRY OPERATOR ELIMINATED)
C
   10 DO 170 I=1,NSYM
        IPOINT(I)=1
  170   I30 = 3
        CALL RECENT(T(1,I),I30)
C
C   PREPARE THE TRANSLATIONAL SYMMETRY MATRIX TU BY ADDING THE
C   APPROPRIATE LATTICE CENTERING TRANSLATIONS
C
      IF(NTT.EQ.1) GO TO 100
      DO 20 II=2,NTT
   20   CALL LTRANS(NLT,TA(1,II),TU)
  100 IF (NDEBUG.NE.0) WRITE (NDEBUG,2100) GLAT,GCFL,NLT,((TU(I1,I2),I2
     1  =1,14),I1=1,3)
 2100 FORMAT (/,1X,'++DEBUG++SPOST++2100++ GLAT = ',A1,', GCFL = ',A1,
     1  ', NLT =',I3,3(/,5X,14F5.2))
C
C   SCAN THROUGH ALL THE SYMMETRY OPERATORS
C
      DO 700 II=2,NSYM
        I=II
        IF(IPOINT(I).NE.1) GOTO 700
        IPOINT(I)=0
C
C   ELIMINATE SCREW AXIS AND GLIDE PLANES WHEN THEY CANNOT GENERATE A SPEC.POS.
C
        IF(GLAT.EQ.'I'.OR.GLAT.EQ.'F'.OR.GLAT.EQ.'H') GOTO 105
        IF(S(1,1,I).GT.0.01.AND.T(1,I).GT.0.01) THEN
          IF(INTEGER(T(1,I),1).EQ.5.AND.
     1      (GLAT.EQ.'B'.OR.GLAT.EQ.'C')) GOTO 190
          GOTO 700
          ENDIF
  190   IF(S(2,2,I).GT.0.01.AND.T(2,I).GT.0.01) THEN
          IF(INTEGER(T(2,I),1).EQ.5.AND.
     1      (GLAT.EQ.'A'.OR.GLAT.EQ.'C')) GOTO 192
          GOTO 700
          ENDIF
  192   IF(S(3,3,I).GT.0.01.AND.T(3,I).GT.0.01) THEN
          IF(INTEGER(T(3,I),1).EQ.5.AND.
     1      (GLAT.EQ.'B'.OR.GLAT.EQ.'A')) GOTO 105
          GOTO 700
          ENDIF
C
C   IF SYMOP OR PRODUCT OF SYMOPS ELIMINATED, DROP IT (SEE STAT. 540)
C
  105   N=II
        IF(I.EQ.II) GOTO 106
        IF(IPOINT(MULT(I,II)).NE.1) GOTO 500
        N=MULT(I,II)
        IPOINT(N)=0
  106   I30 = 3
        CALL MMINM(SD,S(1,1,1),S(1,1,N),I30)
C
C   CHECK FOR A TWO-FOLD AXIS AND SET 222-FLAG TO CREATE THE POSITION
C   AT THE INTERSECTION OF THREE TWO-FILD AXES
C
        I222L=I222
        IF(I222.GT.2) GOTO 107
        DO 250 M1=1,3
        IF(S(M1,M1,N).GT.0.01.AND.T(M1,N).GT.0.01) GOTO 107
        IF(ABS(S(M1,M1,N)).LT.0.99.OR.ABS(S(M1,M1,N)).GT.1.01) GOTO 107
        DO 250 M2=1,3
        IF(M2.EQ.M1) GOTO 250
        IF(ABS(S(M1,M2,N)).GT.0.01) GOTO 107
  250   CONTINUE
        TRACE=S(1,1,N)+S(2,2,N)+S(3,3,N)
        IF(TRACE.GT.-.99.OR.TRACE.LT.-1.01) GOTO 107
        I222=I222+1
        IPOINT(N)=-1
        DO 255 M1=1,3
  255   IF(R222(M1).LT.0.01) R222(M1)=T(M1,N)
        IF(NDEBUG.NE.0) WRITE (NDEBUG,2110) I222
 2110   FORMAT(' ++DEBUG++SPOST++2110++ I222 =',I4)
C
C   SKIP IF MATRIX IS NUL (PURE TRANSLATION)
C
  107   DO 110 K=1,3
        TS(K)=T(K,N)
  110   AROW(K)=ABS(SD(K,1))+ABS(SD(K,2))+ABS(SD(K,3))
        IF((AROW(1)+AROW(2)+AROW(3)).LT.0.01) GOTO 700
C
C   SCAN THROUGH ALL THE LATTICE TRANSLATIONS STORED IN =TU=
C
        NSP3=NSP
        ICUB1=0
  605   DO 600 J=1,NLT
        IF(J.LT.I2LIM) GOTO 600
        IF(J.GT.8.AND.GLAT.EQ.'H') GOTO 600
          CALL VPLUSV(TD,TS,TU(1,J))
C
C   SKIP IF MATRIX DOES NOT LEAD TO SPECIAL POSITION
C
          DO 150 K=1,3
            IF(TD(K).GT.1.99) GOTO 600
            AROW(K)=ABS(SD(K,1))+ABS(SD(K,2))+ABS(SD(K,3))
            IF(AROW(K).LT.0.01) THEN
              IF(TU(K,J).GT.0.99) GOTO 600
              IF(TD(K).GT.0.99) TD(K)=TD(K)-1.0
              IF(TD(K).GT.0.01) GOTO 600
              ENDIF
  150       CONTINUE
          IF(NDEBUG.NE.0) THEN
            WRITE(NDEBUG,2001) SECOND(3),II,I,N,J
 2001       FORMAT(' ++DEBUG++SPOST++2001++(',F6.2,') SYMOP=',I3,',',I3
     +        ,',',I3,', TRANSL=',I3)
            WRITE(NDEBUG,2002)((SD(K,L),L=1,3),TD(K),T(K,N),AROW(K),
     +        K=1,3)
 2002       FORMAT(' ++DEBUG++SPOST++2002++      SD',8X,'TD',8X,'T',
     +        7X,'AROW',3(/,18X,3F5.2,3(4X,F5.2)))
            ENDIF
C
C   TO SOLVE THE EQUATION (SD)(X)=(TD) WE HAVE TO TEACH THE COMPUTER
C   SOME ALGEBRA, BECAUSE SD MAY BE SINGULAR OR NON-QUADRATIC
C
C   SET THREE LINE-MEMORIES
C
          MEM=0
          MEM1=0
          MEM2=0
C
C   INCREMENT NSP AND SET SOME STATISTICS
C
          NSP=NSP+1
          IF (NSP.GT.MSP) GOTO 160
          ICUB=1
          TRACE=0.0
          DO 200 K=1,3
            TP(K,NSP)=0.0
            TRACE=TRACE+ABS(SD(K,K))
            IF(INTEGER(SD(K,K),0).NE.1) ICUB=0
            LR(K)=0
            L2(K)=0
            LC(K)=0
            ROW(K)=0.0
            COL(K)=0.0
            ACOL(K)=0.0
            DO 210 L=1,3
              SP(K,L,NSP)=0.0
              IF(ABS(SD(K,L)).GT.0.01) LR(K)=LR(K)+1
              IF(ABS(SD(L,K)).GT.0.01) LC(K)=LC(K)+1
              IF(ABS(SD(K,L)).GT.1.01) L2(K)=L2(K)+1
              ROW(K)=ROW(K)+SD(K,L)
              COL(K)=COL(K)+SD(L,K)
              ACOL(K)=ACOL(K)+ABS(SD(L,K))
  210         CONTINUE
  200       CONTINUE
          IF(NDEBUG.NE.0) THEN
            WRITE(NDEBUG,2020) (K,LR(K),LC(K),L2(K),ROW(K),AROW(K),
     1        COL(K),ACOL(K),ICUB,K=1,3)
 2020       FORMAT(' ++DEBUG++SPOST++2020   K  LR  LC  L2   ROW  AROW',
     1        '   COL  ACOL  ICUB',3(/,21X,4I4,4F6.2,I4))
            ENDIF
C
C   LOOP THROUGH ALL ROWS...
C
          DO 220 K=1,3
            IF(ACOL(K).LT.0.01.AND.AROW(K).LT.0.01) THEN
C   --        NON-FIXED AXIS
              SP(K,K,NSP)=1.00
              GOTO 220
              ENDIF
            IF(LR(K).EQ.1.AND.LC(K).EQ.1) THEN
C   --        FIXED TRANSLATION
              IF(TD(K).LT.0.01) GOTO 220
              IF(ABS(SD(K,K)).LT.0.01) GOTO 220
              TP(K,NSP)=TD(K)/SD(K,K)
              GOTO 220
              ENDIF
            IF(L2(K).NE.0.AND.ACOL(K).LT.2.01) THEN
C   --        SOMETHING LIKE    2X - Y = 0, T
              SP(K,K,NSP)=1.0
              DO 230 L=1,3
                IF(SD(K,L).LT.-0.01) THEN
                  IF(L.LT.K) THEN
C   --            E.G. X DEPENDENT ON Y, DROP
                    NSP=NSP-1
                    IF(NDEBUG.NE.0) WRITE(NDEBUG,2180)
 2180               FORMAT(/,' ++DEBUG++SPOST++2180++ DEPENDENT AXIS,',
     1                ' OPERATOR DROPPED',/)
                    GOTO 700
                    ENDIF
                  SP(L,K,NSP)=2.0
                  ENDIF
  230           CONTINUE
              IF(TD(K).GT.0.99) TD(K)=TD(K)-1.0
              TP(K,NSP)=TD(K)/2.0
              GOTO 220
              ENDIF
            IF(LR(K).EQ.2.AND.LC(K).EQ.2) THEN
C   --        SOMETHING LIKE    X - Y = 0
              IF(TD(K).GT..99) TD(K)=TD(K)-1.0
              IF(MEM.EQ.0) THEN
C   --          NON-FIXED AXIS OR PURE TRANSLATION ??
                IF(ICUB.EQ.1) GOTO 240
                IF(INTEGER(ROW(K),2).EQ.INTEGER(COL(K),2)) GOTO 240
C   --          PURE TRANSLATION
                IF(MEM1.EQ.0) THEN
C   --            SAVE LINE NUMBER FOR TRANSLATION
                  MEM1=K
                ELSE IF(INTEGER(TD(MEM1),2).EQ.INTEGER(TD(K),2)) THEN
C   --            ZERO OR 1/2 TRANSLATION FOR LINES MEM1 AND K
                  TP(MEM1,NSP)=0.0
                  TP(K,NSP)=0.0
                  IF(COL(K).GT.0.01) TP(K,NSP)=TD(K)
                  IF(COL(MEM1).GT.0.01) TP(MEM1,NSP)=TD(MEM1)
                  IF(I11Z.EQ.0) THEN
                    IF(L2(MEM1).EQ.0.AND.L2(K).EQ.0) THEN
C   --              TETRAGONAL, SET 11Z-FLAG
                      IF(TD(K).LT.0.01) THEN
                        I11Z=1
                        I11ZS=10*MEM1+K
                        ENDIF
                    ELSE IF(INTEGER(TD(MEM1),3).EQ.INTEGER(TD(K),3).
     1                AND.TD(K).GT.0.01) THEN
C   --              GARBAGE
                      NSP=NSP-1
                      IF(NDEBUG.NE.0) WRITE(NDEBUG,2190)
 2190                 FORMAT(/,' ++DEBUG++SPOST++2190++ THIS TRANS',
     1                  'LATION IS NONSENSE, DROPPED',/)
                      GOTO 600
                    ELSE
C   --              HEXAGONAL, E.G. 2X+Y=0
                      TP(MEM1,NSP)=1./3.
                      TP(K,NSP)=1./3.
                      IF(L2(MEM1).NE.0) TP(K,NSP)=2./3.
                      IF(L2(K).NE.0) TP(MEM1,NSP)=2./3.
                      IF(IHEXZ.EQ.0) THEN
                        IHEXZ=1
                        IHEXS=10*MEM1+K
                      ELSE IF(IHEXE.NE.0) THEN
                        TP(MEM1,NSP)=0.0
                        TP(K,NSP)=0.0
                        ENDIF
                      ENDIF
                    ENDIF
                ELSE
C   --            NON-ZERO TRANSLATION FOR LINES MEM1 AND K
                  K1=MEM1
                  KK=K
                  IF(ROW(MEM1).LT.1.01) THEN
                    MEM1=K
                    KK=K1
                    ENDIF
                  TP(MEM1,NSP)=(TD(MEM1)+TD(KK))/2.0
                  TP(KK,NSP)=(TD(MEM1)+TD(KK))/2.0
                  IF(COL(MEM1).LT.0.01) TP(MEM1,NSP)=(TD(MEM1)-TD(KK))
     1              /2.0
                  IF(COL(KK).LT.0.01) TP(KK,NSP)=(TD(MEM1)-TD(KK))/2.0
                  IF(L2(KK).NE.0) THEN
                    TP(KK,NSP)=TD(MEM1)+TD(KK)
                    IF(TP(KK,NSP).GT.0.49.AND.TP(KK,NSP).LT.0.99)
     1                TP(KK,NSP)=TP(KK,NSP)+1.0
                    TP(KK,NSP)=TP(KK,NSP)/3.0
                    TP(MEM1,NSP)=TD(MEM1)-TP(KK,NSP)
                    ENDIF
                  MEM1=K1
                  ENDIF
                GOTO 220
C   --          NON-FIXED AXIS, SAVE LINE NUMBER
  240           MEM=K
                SP(K,K,NSP)=1.0
              ELSE IF(COL(K).GT.1.01.OR.ROW(K).GT.1.01) THEN
C   --          E.G. Y = -X + T
                IF(ICUB.NE.0) THEN
C   --            ONLY FOR NON-CUBICS
                  NSP=NSP-1
                  IF(NDEBUG.NE.0) WRITE(NDEBUG,2250)
 2250             FORMAT(/,' ++DEBUG++SPOST++2250++ DEPENDENT AXIS,',
     1              ' OPERATOR DROPPED',/)
                  IF(I.EQ.II) GOTO 700
                  GOTO 540
                  ENDIF
                I11Z=1
                I11ZS=10*MEM+K
                SP(K,MEM,NSP)=-1.0
                TP(MEM,NSP)=TD(K)
                TP(K,NSP)=ABS(TD(MEM)-TD(K))
                SUMM=ROW(MEM)+ROW(K)+COL(MEM)+COL(K)
                IF(SUMM.GT.7.1.AND.INTEGER(TD(MEM),3).NE.INTEGER(TD(K),
     1            3)) THEN
C   --            X+X .NE. X+X !!
                  I11Z=0
                  I11ZS=0
                  NSP=NSP-1
                  IF(NDEBUG.NE.0) WRITE(NDEBUG,2255)
 2255             FORMAT(/,' ++DEBUG++SPOST++2255++ NO SOLUTION FOR ',
     1              'THIS TRANSLATION',/)
                  GOTO 600
                ELSE IF(INTEGER(TD(MEM),3).EQ.INTEGER(TD(K),3).AND.
     1            INTEGER(TD(K),2).NE.0.AND.INTEGER(TD(K),2).NE.50)
     2            THEN
                  I11Z=0
                  I11ZS=0
                  NSP=NSP-1
                  IF(NDEBUG.NE.0) WRITE(NDEBUG,2255)
                  GOTO 600
                  ENDIF
              ELSE
C   --          E.G. Y = X + T
                I11Z=1
                I11ZS=10*MEM+K
                SP(K,MEM,NSP)=1.0
                TP(K,NSP)=TD(K)
                IF(INTEGER(TD(MEM),3).EQ.INTEGER(TD(K),3).AND.
     1            INTEGER(TD(K),2).NE.0.AND.INTEGER(TD(K),2).NE.50)
     2            THEN
                  I11Z=0
                  I11ZS=0
                  NSP=NSP-1
                  IF(NDEBUG.NE.0) WRITE(NDEBUG,2290)
 2290             FORMAT(/,' ++DEBUG++SPOST++2290++ NO SOLUTION FOR ',
     1              'THIS TRANSLATION',/)
                  GOTO 600
                  ENDIF
                ENDIF
              GOTO 220
              ENDIF
            IF(AROW(K).GT.0.01) THEN
C   --      MAYBE SOMETHING LIKE X=1/3, 2X=2/3
              IF(TD(K).GT.0.99) TD(K)=TD(K)-1
              IF(MEM2.EQ.0) THEN
                MEM2=K
              ELSE IF(L2(MEM2).EQ.0) THEN
C   --        SKIP OPERATOR WITH E.G. X DEPENDENT ON Y
                IF(ABS(SD(MEM2,MEM2)).LT.0.01) THEN
                NSP=NSP-1
                IF(NDEBUG.NE.0) WRITE(NDEBUG,2170)
 2170           FORMAT(/,' ++DEBUG++SPOST++2170++ DEPENDENT AXIS,',
     1            ' OPERATOR DROPPED',/)
                IF(I.EQ.II) GOTO 700
                GOTO 540
                ENDIF
              ELSE IF(L2(K).EQ.0.AND.L2(MEM2).EQ.1) THEN
                IF(TD(K).LT.0.01.OR.INTEGER(TD(MEM2),2).NE.INTEGER
     1            (TD(K),2)) THEN
C   --          AN OTHER E.G. 00Z OR GARBAGE, DROP IT
                  NSP=NSP-1
                  IF(NDEBUG.NE.0) WRITE(NDEBUG,2240)
 2240             FORMAT(/,' ++DEBUG++SPOST++2240++ NON-SPECIFIC, OPE',
     1              'RATOR DROPPED',/)
                  IF(I.EQ.II) GOTO 700
                  GOTO 540
                  ENDIF
                TP(MEM2,NSP)=TD(K)
                ENDIF
              ENDIF
  220    CONTINUE
         IF(I11Z.EQ.1) THEN
C   --   SET FLAG FOR TESTS IN CUBIC CASES
           M1=I11ZS/10
           M2=I11ZS-10*M1
           MM=6-M1-M2
           IF(SD(MM,MM).GT.1.99) I11ZC=MM
           IF(ABS(SD(MM,MM)).LT.0.01) I11ZC=-MM
           ENDIF
         IF(ICUB.EQ.1) THEN
C   --   CUBIC X,X,X, ONLY ZERO-SHIFTS REASONABLE
           TP(1,NSP)=0.0
           TP(2,NSP)=0.0
           TP(3,NSP)=0.0
         ELSE
C   --   SHIFT TP INTO RANGE [0,.99]
           I30 = 3
           CALL RECENT(TP(1,NSP),I30)
           ENDIF
         IF(NDEBUG.NE.0) THEN
           WRITE(NDEBUG,2030) ((SP(K,L,NSP),L=1,3),TP(K,NSP),K=1,3)
 2030      FORMAT(' ++DEBUG++SPOST++2030++ NON-SHIFTED SOLUTION FOUND',
     1       3(/,18X,3F5.2,4X,F6.3))
           ENDIF
         NSP5=NSP-1
         DETERM=0.0
C
C   SHIFT TP CLOSEST TO ORIGIN
C
  265    DETER2=-DETERM
         DO 260 K=1,3
           AROW(K)=ABS(SP(K,1,NSP))+ABS(SP(K,2,NSP))+ABS(SP(K,3,NSP))
           DETERM=DETERM+AROW(K)
           KZZ(K)=INTEGER(AROW(K),0)
  260      IF(KZZ(K).GT.0) KZZ(K)=1
         DETER2=DETER2+DETERM
         CALL REDEFI(SP(1,1,NSP),TP(1,NSP))
         CALL NEWAXS(SP(1,1,NSP),TP(1,NSP))
         IF(NLT.LE.8) GOTO 360
         IF(GLAT.EQ.'H'.AND.SP(3,3,NSP).LT.0.01) GOTO 360
         CALL SHCTOR(TP(1,NSP),SP(1,1,NSP),TU,NLT)
         IF(NDEBUG.NE.0) WRITE(NDEBUG,2035) ((SP(K,L,NSP),L=1,3),
     1     TP(K,NSP),K=1,3)
 2035    FORMAT(' ++DEBUG++SPOST++2035++ SHIFTED SPECIAL POSITION',
     1     3(/,18X,3F5.2,4X,F6.3))
C
C   CHECK IF IT DOES NOT OBVIOUSLY DUPLICATE A PREV. ONE
C
  360    NSP1=NSP-1
         DO 400 K=1,NSP1
           IF(ICOMMA(SP(1,1,NSP),SP(1,1,K),3,2).EQ.0.AND.
     1       ITESTR(TP(1,NSP),TP(1,K),TU,NLT,SP(1,1,K)).EQ.0) THEN
             NSP=NSP1
             LABL=400
             IF(NDEBUG.NE.0) WRITE(NDEBUG,2300) LABL,K
             GOTO 640
             ENDIF
  400     CONTINUE
C
C   CALCULATE SYMMETRY RELATED SPECIAL POSITIONS
C
          ICUB2=0
          DO 440 K=1,NSYM
            ISPTAB(K,NSP)=K
            CALL SYMULT(SPS(1,1,K),TPS(1,K),S(1,1,K),T(1,K),
     +         SP(1,1,NSP),TP(1,NSP))
            I30 = 3
            CALL RECENT(TPS(1,K),I30)
         IF(NDEBUG.NE.0) WRITE(NDEBUG,2350) SECOND(3),K,((SPS(LL,L,K),
     1     L=1,3),TPS(LL,K),LL=1,3)
 2350    FORMAT(' ++DEBUG++SPOST++2350++(',F6.2,') SEITZ MATRIX NR.',I4,
     1     ' OF SPECIAL POSITION',3(/,10X,3F5.2,F8.3))
            IF(I11ZC.GT.0.AND.INTEGER(SPS(I11ZC,I11ZC,K),1).NE.0
     1        .AND.I11ZC.NE.3) THEN
C   --        DEPENDENT TWO-FOLD AXIS IN CUBIC SYSTEM
              NSP=NSP1
              I11ZC=0
              IF(NDEBUG.NE.0) WRITE(NDEBUG,2260)
 2260         FORMAT(/,' ++DEBUG++SPOST++2260++ DEPENDENT AXIS,',
     1          ' OPERATOR DROPPED',/)
              IF(I.EQ.II) GOTO 700
              GOTO 540
            ELSE IF(I11ZC.LT.0) THEN
C   --        DEPENDENT TWO-FOLD AXIS IN CUBIC SYSTEM, TOO
              KK=ABS(I11ZC)
              IF(INTEGER(SPS(KK,KK,K),2).NE.1.AND.KK.NE.3) THEN
                NSP=NSP1
                I11ZC=0
                IF(NDEBUG.NE.0) WRITE(NDEBUG,2260)
                IF(I.EQ.II) GOTO 700
                GOTO 540
                ENDIF
              ENDIF
            DO 375 L=1,3
             KZ(L)=0
             DO 370 LL=1,3
  370          IF(ABS(SPS(L,LL,K)).GT.0.01) KZ(L)=1
C
C   IN CASE E.G. Y IS SYMMETRICALLY EQUIVALENT TO X DROP POSITION
C   WITH Y UNIQUE (PROVES TO BE QUITE TRICKY...)
C
             IF(L.EQ.ICRFL3.AND.ICRFL2.EQ.N) GOTO 375
             IF(KZ(L).EQ.0.AND.KZZ(L).NE.0) THEN
               IF(NDEBUG.NE.0) WRITE(NDEBUG,2150) L,ICRFL1,ICRFL2,
     1           ICRFL3,ICUB2,SPS(1,1,1)
 2150          FORMAT(/,' ++DEBUG++SPOST++2150++ L, ICRFL1,2,3,',
     1           ' ICUB2: ',5I4,', SPS(1,1,1): ',F5.2,/)
               IF(ICRFL2.NE.N) THEN
                 ICRFL2=N
                 IF(ICRFL1.LE.0) THEN
                   ICRFL1=NSP
                   ICRFL3=L
                   IF((KZZ(1)+KZZ(2)+KZZ(3)).EQ.1) THEN
                     DO 475 M1=1,3
                     IF((ABS(SPS(M1,1,K))+ABS(SPS(M1,2,K))+ABS(SPS(M1,3,
     1                 K))).GT.0.01) KZZ(M1)=1
  475                CONTINUE
                     ICUB2=1
                     ENDIF
                 ELSE IF(L.LT.ICRFL3) THEN
                   DO 470 M1=1,3
                     TP(M1,ICRFL1)=TP(M1,NSP)
                     DO 470 M2=1,3
  470                SP(M1,M2,ICRFL1)=SP(M1,M2,NSP)
                   NSP=NSP1
                   ICRFL3=L
                   IF(NDEBUG.NE.0) WRITE(NDEBUG,2130) ICRFL1
 2130              FORMAT(/,' ++DEBUG++SPOST++2130++ SPEC.POS. NR',
     1               I4,' REPLACED BY CURRENT ONE',/)
                   I11Z=0
                   IF(I222.EQ.3) NSP3=0
                   GOTO 500
                 ELSE IF(L.GT.ICRFL3) THEN
                   NSP=NSP1
                   IF(NDEBUG.NE.0) WRITE(NDEBUG,2140) ICRFL1
 2140              FORMAT(/,' ++DEBUG++SPOST++2140++ DEPENDENT AXIS,'
     1               ,' CURRENT SPEC.POS. EQUIVALENT TO NR.',I4,/)
                   I11Z=0
                   IF(I222.EQ.3) NSP3=0
                   GOTO 500
                   ENDIF
               ELSE
                 ICCH=6-ICRFL3-L
                 IF(ICUB2.EQ.1) THEN
                   ICCH=0
                   IF(SPS(1,1,1).GT.0.01) ICCH=4
                   ENDIF
                 IF(ICCH.LT.ICRFL3.OR.ICCH.LT.L) THEN
                   IF(NDEBUG.NE.0) WRITE(NDEBUG,2230)
 2230              FORMAT(/,' ++DEBUG++SPOST++2230++ DEPENDENT AXIS',
     1               ', OPERATOR DROPPED',/)
                   IF(ICRFL1.EQ.NSP) ICRFL1=0
                   NSP=NSP1
                   IF(I11Z.GT.1) THEN
                     I11Z=0
                     GOTO 600
                     ENDIF
                   I11Z=0
                   IF(I222.EQ.3) GOTO 505
                   IF(I.EQ.II) GOTO 590
                   GOTO 500
                   ENDIF
                 IF(ICRFL1.GE.NSP) ICRFL1=-1
                 ENDIF
               ENDIF
  375        CONTINUE
C
C   CALCULATE THE SPECIAL POSITION MULTIPLICATION TABLE =ISPTAB=
C
            K1=K-1
            IF(K1.LE.0) GO TO 440
            DO 420 M1=1,3
C   --        PREVENT REDEFINING OF ORIGINS
              DO 420 M2=1,3
  420         SPP(M1,M2)=0.0
            DO 430 L=1,K1
              IF(ICOMMA(SPS(1,1,L),SPS(1,1,K),3,2).EQ.0.AND.
     +          ITESTR(TPS(1,K),TPS(1,L),TU,NLT,SPP).EQ.0) THEN
                ISPTAB(K,NSP)=L
                GO TO 440
                ENDIF
  430         CONTINUE
  440       CONTINUE
          IF(NDEBUG.NE.0) WRITE(NDEBUG,2440) SECOND(3),NSP,
     +     (ISPTAB(K,NSP),K=1,NSYM)
 2440     FORMAT(/,' ++DEBUG++SPOST++2440++(',F6.2,') NSP=',I3,
     +      ' ISPTAB :',/,10X,48I2)
C
C   GENERATE MULTIPLICITIES OF SPECIAL POSITIONS
C
          ISAVE(1)=ISPTAB(1,NSP)
          L=1
          DO 830 JJ=2,NSYM
            DO 820 K=1,L
  820         IF(ISAVE(K).EQ.ISPTAB(JJ,NSP)) GO TO 830
            L=L+1
            ISAVE(L)=ISPTAB(JJ,NSP)
  830       CONTINUE
          MULTSP(NSP)=L
          IF(L.GE.NSYM) THEN
            NSP=NSP1
            IF(NDEBUG.NE.0) WRITE(NDEBUG,2210)
 2210       FORMAT(/,' ++DEBUG++SPOST++2210++ NSYM:L-CHECK: ABOVE',
     1        ' POSITION NOT SPECIAL, OPERATOR DROPPED',/)
            GOTO 550
            ENDIF
          IF(I11Z.EQ.1.OR.IHEXZ.EQ.1) NSP2=NSP
          IF(I11Z.GT.1) THEN
            IF(MULTSP(NSP2).LE.L) THEN
              NSP=NSP1
              IF(NDEBUG.NE.0) WRITE(NDEBUG,2200)
 2200         FORMAT(/,' ++DEBUG++SPOST++2200++ 11Z-CHECK: ABOVE ',
     1          'SPECIAL POSITION DROPPED',/)
              GOTO 640
              ENDIF
            ENDIF
          IF(I222.NE.I222L) NSP4=NSP
C
C   LOOP THROUGH =ISPTAB=, PICK UP NEW COORDINATES, REDEFINE
C   FREE VARIABLES AND AXES WHERE POSSIBLE
C
          DO 610 K=1,NSYM
            IF(ISPTAB(K,NSP).LT.K) GOTO 610
            ATRAC=0.0
            DO 620 M=1,3
              AROW(M)=ABS(SPS(M,1,K))+ABS(SPS(M,2,K))+ABS(SPS(M,3,K))
              ATRAC=ATRAC+ABS(SPS(M,M,K))
  620         KZ(M)=INTEGER(AROW(M),0)
            IF(K.EQ.1) THEN
              ATRAC2=ATRAC
            ELSE IF(INTEGER(ATRAC,0).EQ.(INTEGER(ATRAC2,0)-2)) THEN
              CALL TEHECU(SP(1,1,NSP),SPS(1,1,K))
              CALL REDEFI(SPS(1,1,K),TPS(1,K))
              CALL NEWAXS(SPS(1,1,K),TPS(1,K))
              IF(NLT.EQ.8) GOTO 630
              IF(GLAT.EQ.'H'.AND.SPS(3,3,K).LT.0.01) GOTO 630
              CALL SHCTOR(TPS(1,K),SPS(1,1,K),TU,NLT)
            ELSE
              CALL REDEFI(SPS(1,1,K),TPS(1,K))
              CALL NEWAXS(SPS(1,1,K),TPS(1,K))
              IF(NLT.EQ.8) GOTO 630
              IF(GLAT.EQ.'H'.AND.SPS(3,3,K).LT.0.01) GOTO 630
              CALL SHCTOR(TPS(1,K),SPS(1,1,K),TU,NLT)
              ENDIF
  630       IF(NDEBUG.NE.0) WRITE(NDEBUG,2280) K,((SPS(I1,I2,K),I2=1,3),
     1        TPS(I1,K),I1=1,3)
 2280       FORMAT(' ++DEBUG++SPOST++2280++ REDEFINED AND SHIFTED ',
     1        'TRIPPLET NR:',I4,3(/,10X,3F5.2,F8.3))
C
C   REMOVE THE CURRENT SPECIAL POSITION IF IT OR ANY SYMMETRY
C   DERIVATIVE OF IT DUPLICATES AN EXISTING SPECIAL POSITION
C
            IF(NSP1.LE.0) GO TO 640
  395       IIFL=0
            DO 390 L=1,NSP1
              IF(MULTSP(L).NE.MULTSP(NSP)) GOTO 390
            DO 480 MM=1,2
C   --        CHECK SPS...
              SSIG=1.0
C   --        ...AS WELL AS -SPS
              IF(MM.EQ.2) SSIG=-1.0
              AOTR=0.0
              DO 490 M1=1,3
                DO 490 M2=1,3
                IF(M1.NE.M2) AOTR=AOTR+ABS(SPS(M1,M2,K))
  490           SPP(M1,M2)=SPS(M1,M2,K)*SSIG
              IF(ICOMMA(SPP,SP(1,1,L),3,2).EQ.0) THEN
  495           IF(ITESTR(TP(1,L),TPS(1,K),TU,NLT,SPS(1,1,K)).EQ.0) THEN
                  NSP=NSP1
                  LABL=495
                  IF(NDEBUG.NE.0) WRITE(NDEBUG,2300) LABL,L
 2300             FORMAT(/,' ++DEBUG++SPOST++',I4,'++ ABOVE POSITION',
     1              ' EQUIVALENT TO NR.',I3,', DROPPED',/)
                  GOTO 640
                ELSE IF(I11Z.EQ.1.AND.ICUB.EQ.0.AND.IIFL.NE.1) THEN
C   --            E.G. (T(X),T(Y),Z): CHECK ALSO (T(Y),T(X),Z)...
C                 AND (0,T(Y)-T(X),Z)
                  M1=1
                  IF(KZ(1).NE.1.OR.SPS(1,1,K).GT.1.01) M1=2
                  M2=M1+1
                  IF(KZ(M2).NE.1.OR.SPS(2,2,K).GT.1.01) M2=M2+1
                  IF(M2.EQ.4.OR.KZ(M2).NE.1.OR.SPS(M2,M2,K).GT.1.01)
     1              THEN
                    IIFL=1
                    GOTO 495
                    ENDIF
                  M=6-M1-M2
                  DO 570 MX=1,2
                  TT(M)=TPS(M,K)
                  TT(M1)=TPS(M2,K)
                  TT(M2)=TPS(M1,K)
                  IF(MX.EQ.2) THEN
                    TT(M1)=0
                    TT(M2)=TPS(M2,K)-TPS(M1,K)
                    I30 = 3
                    CALL RECENT(TT,I30)
                    ENDIF
                  IF(ITESTR(TP(1,L),TT,TU,NLT,SPS(1,1,K)).EQ.0) THEN
                    NSP=NSP1
                    ICRFL1=0
                    LABL=570
                    IF(NDEBUG.NE.0) WRITE(NDEBUG,2300) LABL,L
                    GOTO 640
                    ENDIF
  570             CONTINUE
                  ENDIF
                ENDIF
              IF(ATRAC.LT.0.01.OR.AOTR.LT.0.01) GOTO 455
C   --        NON-DIAGONAL ELEMENTS IN SPS: CHECK ALSO TRANSPOSE OF SPS
              DO 450 M=1,3
              DO 460 M1=1,3
                TPP(M1)=TP(M1,L)
                TSS(M1)=TPS(M1,K)
                DO 460 M2=1,3
                SS(M1,M2)=SPS(M2,M1,K)*SSIG
  460           SPP(M1,M2)=SP(M1,M2,L)
              IF(M.EQ.2) THEN
C   --          A "CUBIC-SPECIAL"
                CALL SYMULT(SPP,TPP,S(1,1,K),T(1,K),SP(1,1,L),TP(1,L))
                I30 = 3
                CALL RECENT(TPP,I30)
              ELSE IF(INTEGER(SPS(1,1,1),0).NE.0.AND.INTEGER(SPS(1,1,K)
     1          ,0).EQ.0) THEN
C   --          AN OTHER "CUBIC SPECIAL"
                DO 405 M1=1,3
                  IF(INTEGER(SPS(1,M1,K),0).EQ.0) GOTO 405
                  IF(INTEGER(SPS(M1,1,K),0).EQ.0) GOTO 405
                  IF(INTEGER(SPS(M1,M1,K),0).NE.0) GOTO 405
                  SS(1,1)=SS(1,M1)
                  SS(1,M1)=0.0
                  SS(M1,M1)=SS(M1,1)
                  SS(M1,1)=0.0
                  TSS(1)=0.0
                  TSS(M1)=0.0
                  GOTO 415
  405             CONTINUE
                GOTO 450
              ELSE
                GOTO 450
                ENDIF
  415         IF(ICOMMA(SS,SPP,3,2).EQ.0.AND.
     +           ITESTR(TPP,TSS,TU,NLT,SS).EQ.0) THEN
                NSP=NSP1
                LABL=415
                IF(NDEBUG.NE.0) WRITE(NDEBUG,2300) LABL,L
                GOTO 640
                ENDIF
  450         CONTINUE
              GOTO 480
  480       CONTINUE
C
C   REMOVE SITES WITH FIXED COORDINATES IF THEY HAVE THE SAME MULTIPLI-
C   CITY AS SITES WITH FREE COORDINATES IN SAME DIRECTION
C
  455       ATRAC1=0
            DETER1=0.0
            DO 380 M=1,3
              ATRAC1=ATRAC1+ABS(SP(M,M,L))
              TPP(M)=TPS(M,K)
              AROW(M)=0.0
              DO 380 MM=1,3
                SPP(M,MM)=SP(M,MM,L)
                DETER1=DETER1+ABS(SPP(M,MM))
                AROW(M)=AROW(M)+ABS(SPP(M,MM))
  380         CONTINUE
            IF(INTEGER(ATRAC1,0).EQ.INTEGER(ATRAC2,0)) GOTO 390
            DO 465 M=1,3
              IF(ATRAC1.GT.(ATRAC2+0.1).AND.ABS(SP(M,M,L)).LT.(ABS(
     1          SPS(M,M,K))-0.1)) GOTO 390
              IF(ATRAC1.LT.(ATRAC2-0.1).AND.ABS(SP(M,M,L)).GT.(ABS(
     1          SPS(M,M,K))+0.1)) GOTO 390
              IF(ABS(SP(M,M,L)).GT.0.9.AND.ABS(SPS(M,M,K)).GT.0.9
     1          .AND.INTEGER(SP(M,M,L),0).NE.INTEGER(SPS(M,M,K),0))
     2          GOTO 390
              DO 465 MM=1,3
                IF(MM.EQ.M) GOTO 465
                IF(INTEGER(SP(M,M,L),0).EQ.INTEGER(SPS(M,M,K),0).AND.
     1            INTEGER(SP(MM,M,L),0).NE.INTEGER(SPS(MM,M,K),0))
     2            THEN
                  KKK=6-MM-M
                  IF(ABS(SP(MM,M,L)).GT.0.9) THEN
                    IF(ABS(SPS(MM,MM,K)).LT.0.1) GOTO 390
                    IF(ABS(SPS(KKK,M,K)).GT.0.9) GOTO 390
                  ELSE IF(ABS(SPS(MM,M,K)).GT.0.9) THEN
                    IF(ABS(SP(MM,MM,K)).LT.0.1) GOTO 390
                    IF(ABS(SP(KKK,M,K)).GT.0.9) GOTO 390
                    ENDIF
                  ENDIF
  465         CONTINUE
            IF(ATRAC1.LT.(ATRAC2-0.1)) GOTO 905
            CALL REDEFI(SPP,TPP)
            IF(NLT.EQ.8) GOTO 385
            IF(GLAT.EQ.'H'.AND.SPP(3,3).LT.0.01) GOTO 385
            CALL SHCTOR(TPP,SPP,TU,NLT)
  385       IF(ICOMVE(TPP,TP(1,L),3,2).EQ.0) THEN
              NSP=NSP1
              IF(NDEBUG.NE.0) WRITE(NDEBUG,2120) L
 2120         FORMAT(/,' ++DEBUG++SPOST++2120++ ABOVE POSITION DROPPED',
     1          ' IN FIXED COORD. TEST AGAINST POS.',I3,/)
              GOTO 640
              ENDIF
  905       DO 900 M=1,3
              TPP(M)=TP(M,L)
              LR(M)=0
              AROW(M)=0.0
              DO 900 MM=1,3
                SPP(M,MM)=SPS(M,MM,K)
                AROW(M)=AROW(M)+ABS(SPS(M,MM,K))
                LR(M)=INTEGER(AROW(M),0)
  900         CONTINUE
            CALL REDEFI(SPP,TPP)
            IF(NLT.EQ.8) GOTO 910
            IF(GLAT.EQ.'H'.AND.SPP(3,3).LT.0.01) GOTO 910
            CALL SHCTOR(TPP,SPP,TU,NLT)
  910       IF(ICOMVE(TPP,TPS(1,K),3,2).EQ.0) THEN
              NSP=NSP-1
              NSP1=NSP1-1
              NSP3=NSP3-1
              IF(NSP2.GT.L) NSP2=NSP2-1
              IF(NSP4.GT.L) NSP4=NSP4-1
              IF(NSP5.GT.L) NSP5=NSP5-1
              IF(NSP6.GT.L) NSP6=NSP6-1
              IF(ICRFL1.GT.L) ICRFL1=ICRFL1-1
              IF(NSP2.EQ.L) NSP2=NSP
              IF(NSP4.EQ.L) NSP4=NSP
              IF(NSP5.EQ.L) NSP5=NSP
              IF(NSP6.EQ.L) NSP6=NSP
              IF(ICRFL1.EQ.L) ICRFL1=NSP
              DO 940 M1=L,NSP
                DO 920 M=1,3
                  TP(M,M1)=TP(M,M1+1)
                  DO 920 MM=1,3
  920             SP(M,MM,M1)=SP(M,MM,M1+1)
                DO 930 M=1,NSYM
  930             ISPTAB(M,M1)=ISPTAB(M,M1+1)
                MULTSP(M1)=MULTSP(M1+1)
  940           CONTINUE
              IF(NDEBUG.NE.0) WRITE(NDEBUG,2160) L
 2160         FORMAT(/,' ++DEBUG++SPOST++2160++ FIXED COORD. TEST: ',
     1          'POS.',I3,' DROPPED, REST DOWN SHIFTED',/)
              GOTO 395
              ENDIF
  390       CONTINUE
  610     CONTINUE
C
C   IN CASE OF TWO FOLD AXIS RUN ALSO WITH ALL THREE COORDINATES FIXED
C
  640     IF(NDEBUG.NE.0) WRITE(NDEBUG,2270) I222,I222L,NSP,NSP1,NSP5,
     1      NSP6,NSP2,I11Z,DETERM
 2270     FORMAT(' ++DEBUG++SPOST++2270++ I222,I222L,  NSP, NSP1,',
     1      ' NSP5, NSP6, NSP2, I11Z, DETERM',/,22X,8I6,F6.2)
          IF(I11Z.NE.0) GOTO 515
          IF(INTEGER(DETERM,0).NE.1.OR.I222.GT.3) GOTO 600
          IF(NSP.EQ.NSP5) THEN
            I222S=0
            GOTO 600
            ENDIF
          IF(I222S.EQ.0) THEN
            DO 270 M=1,3
              LC(M)=1
  270         IF(ABS(SP(M,M,NSP)).LT.0.01) LC(M)=0
            NSP6=NSP
            I222S=1
            STEPER=0.0
            ENDIF
          NSP=NSP+1
C
C   STEP THROUGH THE ORIGINALLY FREE COORDINATE IN STEPS OF .25
C
          DO 280 M=1,3
            TP(M,NSP)=TP(M,NSP6)
            IF(LC(M).EQ.1) TP(M,NSP)=TP(M,NSP6)+.25*STEPER
            IF(TP(M,NSP).GT.0.9) THEN
              NSP=NSP-1
              I222S=0
              GOTO 600
              ENDIF
            DO 290 MM=1,3
  290         SP(M,MM,NSP)=SP(M,MM,NSP6)
            SP(M,M,NSP)=0.0
  280       CONTINUE
          STEPER=STEPER+1.0
          GOTO 265
C
C   THREE-FOLD, FOUR-FOLD OR SIX-FOLD AXIS, THEREFORE
C   RUN ALSO X=Y=0, X=Y=1/8, X=Y=1/4 ...
C   IN ROMBOHEDRAL OR CUBIC CASES USE X=Y=Z
C
  515     IF(I11Z.EQ.6.OR.NSP.EQ.NSP5) THEN
            I11Z=0
            GOTO 600
            ENDIF
          M1=I11ZS/10
          M2=I11ZS-10*M1
          M4=6-M1-M2
          I11=I11Z-1
          NSP=NSP+1
          TP(M1,NSP)=TP(M1,NSP2)+I11*SP(M1,1,NSP2)*.125
          TP(M2,NSP)=TP(M2,NSP2)+I11*SP(M2,1,NSP2)*.125
          TP(M4,NSP)=TP(M4,NSP2)
          DO 530 MM=1,3
            DO 530 MN=1,3
  530       SP(MM,MN,NSP)=SP(MM,MN,NSP2)
          SP(M1,1,NSP)=0.0
          SP(M2,1,NSP)=0.0
          IF(ICUB.EQ.1.OR.ICUB1.EQ.1) THEN
            ICUB1=1
            TP(M4,NSP)=TP(M4,NSP2)+I11*SP(M4,1,NSP2)*.125
            SP(M4,1,NSP)=0.0
            ENDIF
          I11Z=I11Z+1
          I30 = 3
          CALL RECENT(TP(1,NSP),I30)
          GOTO 265
  600     CONTINUE
        I222L=I222
C
C   CHECK SOME SPECIAL CASES
C
  500   IF(NDEBUG.NE.0) WRITE(NDEBUG,2220) I222,I2LIM,IHEXZ,NSP3,NSP2
 2220   FORMAT(/,' ++DEBUG++SPOST++2220++ FLAGS:   I222,I2LIM,',
     1    'IHEXZ, NSP3, NSP2',/,31X,5I6,/)
        IF(NSP.EQ.NSP3) GOTO 550
        IF(I222.EQ.10) GOTO 590
        IF(I222.EQ.4) I222=5
        IF(ATRAC2.GT.0.01.OR.I444.NE.0) GOTO 505
C
C   IN VERY FEW CASES (.25,.25,.75) IS STILL NOT FOUND...
C
          DO 580 M=1,3
  580       IF(INTEGER(TP(M,NSP),2).NE.25.AND.INTEGER(TP(M,NSP),2)
     1        .NE.75) GOTO 505
          DO 585 M=1,3
  585       TS(M)=2.0*TP(M,NSP)
          I444=1
          GOTO 605
  505   IF(I222.EQ.3) THEN
C
C   THREE TWO-FOLD AXES PRESENT, THEREFORE ARTIFICIALLY
C   ADD POINT WITH 222-SYMMETRY (ELSE NOT FOUND E.G. IN P 222)
C
          I222=4
          IF(NSP.EQ.NSP3.AND.I.EQ.II) I222=10
          DO 510 M1=1,3
            DO 520 M2=1,3
  520       SD(M1,M2)=0.0
            TS(M1)=R222(M1)
  510       SD(M1,M1)=2.0
          GOTO 605
          ENDIF
        IF(IHEXZ.NE.0) THEN
C
C   OPERATOR ALONG HEXAGONAL AXIS, THEREFORE INTERCHANGE X AND Y,
C   THEN SET X = Y = 0
C
          M1=IHEXS/10
          M2=IHEXS-10*M1
          M3=6-M1-M2
          IF(IHEXZ.EQ.1) THEN
            IHEXZ=2
            DUMP1=SD(M1,M1)
            SD(M1,M1)=SD(M2,M2)
            SD(M2,M2)=DUMP1
            DUMP1=SD(M1,M2)
            SD(M1,M2)=SD(M2,M1)
            SD(M2,M1)=DUMP1
            GOTO 605
          ELSE IF(IHEXE.EQ.0) THEN
              IHEXE=1
              GOTO 605
          ELSE
            IHEXZ=0
            IHEXE=0
            ENDIF
          ENDIF
  550   IHEXZ=0
        M3=0
  545   IF(I222.EQ.10) THEN
C   --    A "CUBIC SPECIAL"
          I222=4
          I2LIM=8
          DO 555 MM=1,3
  555     TS(MM)=2.0*TP(MM,NSP4)
          GOTO 605
          ENDIF
        I2LIM=1
        ICUB1=0
C
C   CREATE INTERSECTION OF TWO SYMOPS, USE THE MULTIPLICATION TABLE
C
  540   IF(I.GT.2) THEN
          I=I-1
          GOTO 105
        ELSE
          DO 560 J=I,II
  560     IF(IPOINT(MULT(II,J)).EQ.1) IPOINT(MULT(II,J))=0
          ENDIF
        GOTO 700
  590   IPOINT(N)=-1
        IF(I222.EQ.10) GOTO 545
  700   CONTINUE
      NSP=NSP+1
  160 IF(NSP.GT.MSP) THEN
        WRITE(NOUT,2010) MSP
 2010   FORMAT(/,' **ERROR**SPOST**MORE THAN',I4,'SPECIAL POSITIO',
     1    'NS, SEARCH ABORTED',/)
        NSP=0
        IER=1
        GOTO 800
        ENDIF
C
C   FINALLY ADD THE GENERAL POSITION
C
      DO 710 I=1,3
      TP(I,NSP)=0.0
      DO 710 J=1,3
  710 SP(I,J,NSP)=0.0
      DO 720 I=1,3
  720 SP(I,I,NSP)=1.0
      DO 730 I=1,NSYM
  730 ISPTAB(I,NSP)=I
      MULTSP(NSP)=NSYM
      IF (GLAT.EQ.'P') GOTO 800
      MUL=2
      IF (GLAT.EQ.'F') MUL=4
      IF (GLAT.EQ.'H') MUL=3
      DO 740 I=1,NSP
  740 MULTSP(I)=MULTSP(I)*MUL
C
C   RESET DEBUG-SWITCH
C
  800 NDEBUG=NDEB
      RETURN
      END
C
C
C
C
      SUBROUTINE REDEFI(SPS,TPS)
C
C     CALLS RECENT
C
C     REDEFINES INDEPENDENT VARIABLES IN MATRIX =SPS= SO THAT THE
C     ACCORDING TRANSLATIONS IN =TPS= BECOME ZERO
C
C     WRITTEN BY D. ALTERMATT
C     MCMASTER UNIVERSITY, DEC. 1984
C
      DIMENSION SPS(3,3),TPS(3),AROW(3)
C
      DO 100 I=1,3
        AROW(I)=ABS(SPS(I,1))+ABS(SPS(I,2))+ABS(SPS(I,3))
        IF(AROW(I).LT.0.01) GOTO 100
        K1=I-1
        K2=I+1
        IF(K1.LT.1) K1=3
        IF(K2.GT.3) K2=1
        IF(ABS(SPS(I,I)).GT.0.01.AND.(ABS(SPS(I,K1))+ABS(SPS(I,K2)))
     1    .LT.0.01) THEN
          TPS(K1)=TPS(K1)-SPS(K1,I)*TPS(I)/SPS(I,I)
          TPS(K2)=TPS(K2)-SPS(K2,I)*TPS(I)/SPS(I,I)
          TPS(I)=0.0
          IF(SPS(I,I).LT.-.01) THEN
            SPS(I,I)=ABS(SPS(I,I))
            SPS(K1,I)=-SPS(K1,I)
            SPS(K2,I)=-SPS(K2,I)
            ENDIF
          ENDIF
  100   CONTINUE
      I30 = 3
      CALL RECENT(TPS,I30)
      RETURN
      END
C
C
C
C
      SUBROUTINE NEWAXS(SPS,TPS)
C
C     CALLS INTEGER
C
C     REDEFINES THE AXIS FOR INDEPENDENT VARIABLES SO THAT E.G.
C
C           ( 0 1 0)  (TX)                ( 1 0 0)  ( 0)
C           ( 0 0 0)  (TY)    BECOMES     ( 0 0 0)  (TY)
C           (-1 0 0)  (TZ)                ( 0 0 1)  ( 0)
C
C     WRITTEN BY D. ALTERMATT
C     MCMASTER UNIVERSITY, DEC. 1984
C
      DIMENSION SPS(3,3),TPS(3),NEWFL(3)
C
      DO 5 I=1,3
    5   NEWFL(I)=0
C
      DO 10 I=1,3
        XK=ABS(SPS(I,1))+ABS(SPS(I,2))+ABS(SPS(I,3))
        IK=INTEGER(XK,0)
        IF(IK.NE.1.OR.ABS(SPS(I,I)).GT.0.1) GOTO 10
          K1=I-1
          K2=I+1
          IF(K1.LT.1) K1=3
          IF(K2.GT.3) K2=1
          IF(ABS(SPS(I,K1)).GT.0.1.AND.ABS(SPS(K2,K1)).LT.0.1) THEN
            IF(NEWFL(K1).EQ.0.AND.ABS(SPS(K1,K1)).GT.0.1) GOTO 10
            NEWFL(I)=1
            SPS(I,I)=ABS(SPS(I,K1))
            SPS(I,K1)=0.0
            TPS(I)=0.0
            GOTO 10
          ELSE IF(ABS(SPS(I,K2)).GT.0.1.AND.ABS(SPS(K1,K2)).LT.0.1) THEN
            IF(NEWFL(K2).EQ.0.AND.ABS(SPS(K2,K2)).GT.0.1) GOTO 10
            NEWFL(I)=1
            SPS(I,I)=ABS(SPS(I,K2))
            SPS(I,K2)=0.0
            TPS(I)=0.0
            ENDIF
   10   CONTINUE
      RETURN
      END
C
C
C
C
      SUBROUTINE TEHECU(SP,SPS)
C
C     CALLS INTEGER
C
C     REDEFINES A TRANSFORMED SPECIAL POSITION MATRIX =SPS= IN CASES
C     WITH TWO EQUIVALENT AXES, ON THE BASIS OF NON-ZERO DIAGONAL
C     ELEMENTS IN =SP= :
C
C                    ( 0-1 0)           (-1 0 0)            ( 1 0 0)
C     E.G.   =SPS= : ( 1 0 0)  BECOMES  ( 0 1 0)  IF =SP= : ( 0 1 0)
C                    ( 0 0 Z)           ( 0 0 Z)            ( 0 0 Z)
C
C     WRITTEN BY D. ALTERMATT
C     MCMASTER UNIVERSITY, DEC. 1984
C
      DIMENSION SP(3,3),SPS(3,3),KS(3),KP(3)
C
      DO 10 I=1,3
        KS(I)=INTEGER(SP(I,I),0)
   10   KP(I)=INTEGER(SPS(I,I),0)
C
      DO 20 I=1,3
        IF(KS(I).NE.0.AND.KP(I).EQ.0) THEN
          I1=I+1
          DO 30 J=I1,3
            IF(ABS(SPS(I,J)).GT.0.1.AND.ABS(SPS(J,I)).GT.0.1) THEN
              SPS(I,I)=SPS(I,J)
              SPS(I,J)=0.0
              SPS(J,J)=SPS(J,I)
              SPS(J,I)=0.0
              ENDIF
   30       CONTINUE
          ENDIF
   20   CONTINUE
      RETURN
      END
C
C
C
C
      SUBROUTINE SHCTOR(TP,SP,TU,NLT)
C
C     CALLS RECENT, REDEFI
C
C     SHIFT A GIVEN LATTICE VECTOR TP CLOSEST TO THE ORIGIN WITH
C     RESPECT TO THE LATTICE TRANSLATIONS IN TU AND NON-ZERO ELEMENTS
C     IN MATRIX SP
C
C     WRITTEN BY D. ALTERMATT
C     MCMASTER UNIVERSITY, OCTOBER 1984
C
      DIMENSION TP(3),SP(3,3),T1(3),T2(3),TU(3,NLT)
C
      DO 10 I=1,3
   10 T2(I)=TP(I)
      DO 20 I=9,NLT
        DO 30 K=1,3
   30     T1(K)=TU(K,I)+TP(K)
        I30 = 3
        CALL RECENT(T1,I30)
        CALL REDEFI(SP,T1)
        SUMP=TP(1)+TP(2)+TP(3)
        SUM1=T1(1)+T1(2)+T1(3)
        SUM2=T2(1)+T2(2)+T2(3)
        IF(SUM1.LT.(SUMP-0.01).AND.SUM1.LT.(SUM2-0.01)) GOTO 25
        IF(SUMP.LT.(SUM1-0.01).OR.SUM2.LT.(SUM1-0.01)) GOTO 20
        SUMP=SUMP-TP(3)
        SUM1=SUM1-T1(3)
        SUM2=SUM2-T2(3)
        IF(SUM1.LT.(SUMP-0.01).AND.SUM1.LT.(SUM2-0.01)) GOTO 25
        IF(SUMP.LT.(SUM1-0.01).OR.SUM2.LT.(SUM1-0.01)) GOTO 20
        IF(TP(1).LT.(T1(1)-0.01).OR.T2(1).LT.(T1(1)-0.01)) GOTO 20
   25   DO 40 L=1,3
   40     T2(L)=T1(L)
   20   CONTINUE
      DO 50 I=1,3
   50   TP(I)=T2(I)
      RETURN
      END
C
C
C
C
      SUBROUTINE SPMULT
C     -----------------
C
C     WRITTEN BY D. ALTERMATT
C     MCMASTER UNIVERSITY, NOVEMBER 1984
C
C     CALLS ICOMMA, ITESTR, LTRANS, SYMULT, RECENT
C
C     THIS ROUTINE CALCULATES THE MULTIPLICATION TABLE =ISPTAB=
C *** WHICH IS DEPENDENT ON THE CHOICE OF THE SPECIAL POSITION REPRESENTATIVE.
C *** =ISPTAB= THEREFORE MUST BE RECALCULATED WHENEVER ANY SPECIAL POSITION
C *** REPRESENTATIVE IS REDEFINED !!
C
      IMPLICIT CHARACTER*8 (G)
      IMPLICIT CHARACTER*80 (H)
C
C     SYSTEM COMMON BLOCKS
C
      COMMON/FILES/NIN,NOUT,NDEBUG
      COMMON/SYM/NSYM,MSYM,S(3,3,48),T(3,48),NTT,TA(3,4)
      COMMON/GSYM/GLAT,GCENT,HMSYM,HALL,GSYSP(2,30)
      COMMON/SYMA/MULT(48,48),IT(3,48),NSP,MSP,ISPTAB(48,30),
     +   MULTSP(30),SP(3,3,30),TP(3,30)
C
      DIMENSION TU(3,14),SPS(3,3,48),TPS(3,48),KZ(3),ISAVE(48),SPP(3,3)
C   GENERAL LATTICE TRANSLATIONS
      DATA TU/1.,1.,1., 1.,0.,0., 0.,1.,0., 0.,0.,1., 0.,1.,1.,
     + 1.,0.,1., 1.,1.,0., 0.,0.,0., 18*.0/
C
C   PREPARE THE TRANSLATIONAL SYMMETRY MATRIX TU BY ADDING THE
C   APPROPRIATE LATTICE CENTERING TRANSLATIONS
C
      NLT=8
      IF(NTT.EQ.1) GOTO 10
      DO 5 I=2,NTT
    5   CALL LTRANS(NLT,TA(1,I),TU)
C
C   SET MULTIPLICATION OF POSITION DUE TO LATTICE CENTRINGS
C
   10 MUL=1
      IF(GLAT.EQ.'P') GOTO 50
      MUL=2
      IF(GLAT.EQ.'F') MUL=4
      IF(GLAT.EQ.'H') MUL=3
   50 CONTINUE
C
C   TO CHECK WHETHER THE GENERAL POSITION IS PRESENT, SET A FLAG
C
      IGENFL=0
C
C   MULTIPLY ALL SPECIAL POSITIONS...
C
      DO 100 I=1,NSP
C
C   ...BY ALL SYMMETRY OPERATORS
C
         DO 150 J=1,NSYM
C
C   SET =ISPTAB= ELEMENT
C
            ISPTAB(J,I)=J
            CALL SYMULT(SPS(1,1,J),TPS(1,J),S(1,1,J),T(1,J),SP(1,1,I),
     1           TP(1,I))
            I30 = 3
            CALL RECENT(TPS(1,J),I30)
            IF(J.EQ.1) GOTO 150
C
C   CHECK IF A PREVIOUS MULTIPLICATION GAVE THE SAME RESULT
C
            J1=J-1
            DO 180 K=1,3
C   --        PREVENT REDEFINING OF ORIGINS
              DO 180 K1=1,3
  180         SPP(K,K1)=0.0
            DO 200 K=1,J1
               IF(ICOMMA(SPS(1,1,K),SPS(1,1,J),3,2).EQ.0.AND.
     1            ITESTR(TPS(1,J),TPS(1,K),TU,NLT,SPP).EQ.0) THEN
                  ISPTAB(J,I)=K
                  GOTO 150
                  ENDIF
  200          CONTINUE
  150       CONTINUE
C
C   NOW CALCULATE THE MULTIPLICITY FOR THIS POSITION
C
         ISAVE(1)=ISPTAB(1,I)
         L=1
         DO 250 J=2,NSYM
            DO 300 K=1,L
  300       IF(ISAVE(K).EQ.ISPTAB(J,I)) GOTO 250
            L=L+1
            ISAVE(L)=ISPTAB(J,I)
  250       CONTINUE
         IF(L.EQ.NSYM) IGENFL=1
         IF(NDEBUG.NE.0) WRITE(NDEBUG,2000) GLAT,MUL,NSYM,L
 2000    FORMAT(/,' ++DEBUG++SPMULT++2000++ GLAT = ',A1,', MUL =',
     1     I3,', NSYM =',I3,', L =',I3)
         MULTSP(I)=L*MUL
  100    CONTINUE
C
C   ADD THE GENERAL POSITION IF IGENFL.EQ.0
C
      IF(IGENFL.EQ.0) THEN
         NSP=NSP+1
         DO 320 I=NSP,2,-1
            DO 330 J=1,3
              TP(J,I)=TP(J,I-1)
              DO 330 K=1,3
  330         SP(J,K,I)=SP(J,K,I-1)
            GSYSP(1,I)=GSYSP(1,I-1)
            GSYSP(2,I)=GSYSP(2,I-1)
            MULTSP(I)=MULTSP(I-1)
            DO 340 J=1,NSYM
 340          ISPTAB(J,I)=ISPTAB(J,I-1)
 320        CONTINUE
         DO 350 I=1,3
            TP(I,1)=0.0
            DO 360 J=1,3
  360       SP(I,J,1)=0.0
  350       SP(I,I,1)=1.0
         MULTSP(1)=MUL*NSYM
         DO 400 I=1,NSYM
  400    ISPTAB(I,1)=I
         ENDIF
      RETURN
      END
C
C
C
C
      SUBROUTINE SPESOR
C     -----------------
C
C     SORTS THE SPECIAL POSITIONS FIRST INTERNALLY AND THEN
C     ACCORDING TO THE FOLLOWING RULES:
C
C     A)  BY DECREASING MULTIPLICITY
C     B)  BY NON-FIXED PARAMETERS X,Y,Z
C     C)  BY MINIMUM TRANSLATION IN X,Y,Z-DIRECTION
C
C *** USES A PRELIMINARY VERSION OF =ISPTAB= TO SAVE CALCULATION TIME
C
C     CALLS SYMULT, RECENT, INTEGER, ICOMMA, SHCTOR, REDEFI,
C           TEHECU, NEWAXS
C
C***** CALLS TIME-FUNCTION =SECOND(I)=
C
C     WRITTEN BY DANIEL ALTERMATT
C     MCMASTER UNIVERSITY, OCTOBER 1984
C
      IMPLICIT CHARACTER*8 (G)
      IMPLICIT CHARACTER*80 (H)
      DIMENSION TU(3,14),SPS(3,3,48),TPS(3,48),TA(3),IPT(48),ISS(48)
      DIMENSION ROW(3), GTEMP(2), SS(3,3)
C
C     SYSTEM COMMON BLOCKS
C
      COMMON /FILES/NIN,NOUT,NDEBUG
      COMMON /SYMA/ MULT(48,48),IT(3,48),NSP,MSP,ISPTAB(48,30),
     1   MULTSP(30),SP(3,3,30),TP(3,30)
      COMMON /SYM/  NSYM,MSYM,S(3,3,48),T(3,48),NTT,TNT(3,4)
      COMMON /GSYM/ GLAT,GCENT,HMSYM,HALL,GSYSP(2,30)
C   GENERAL LATTICE TRANSLATIONS
      DATA TU/1.,1.,1., 1.,0.,0., 0.,1.,0., 0.,0.,1., 0.,1.,1.,
     + 1.,0.,1., 1.,1.,0., 0.,0.,0., 18*.0/
C
      NLT=8
      IF(NTT.EQ.1) GOTO 8
      DO 5 I=2,NNT
    5   CALL LTRANS(NLT,TNT(1,I),TU)
C
    8 IF(NDEBUG.NE.0) WRITE(NDEBUG,2100) SECOND(3)
 2100 FORMAT(//,' ++DEBUG++SPESOR++2100++ START AT',F8.3,/)
C
C     FOR ALL SPECIAL POSITIONS DO...
C
      DO 10 N=1,NSP
      STIME=SECOND(3)
C
C     FIRST EXPAND THE SPECIAL POSITIONS TO THEIR FULL SET
C     OF SEITZ MATRICES...
C
      ATRAC1=0.0
      DO 20 I=1,3
        TPS(I,1)=TP(I,N)
        ATRAC1=ATRAC1+ABS(SP(I,I,N))
        ROW(I)=ABS(SP(I,1,N))+ABS(SP(I,2,N))+ABS(SP(I,3,N))
        DO 20 J=1,3
   20     SPS(I,J,1)=SP(I,J,N)
      IPS=1
      IPT(1)=1
C
C     ASSUME GIVEN POSITION TO BE LEGAL, CONVERT IT TO OUR STANDARD
C
      CALL REDEFI(SPS(1,1,1),TPS(1,1))
      CALL NEWAXS(SPS(1,1,1),TPS(1,1))
      IF (NLT.EQ.8) GOTO 22
      IF (GLAT.EQ.'H'.AND.SPS(3,3,1).LT.0.01) GOTO 22
      CALL SHCTOR(TPS(1,1),SPS(1,1,1),TU,NLT)
C
C     NOW EXPAND SET
C
   22 DO 30 I=1,NSYM
      DO 40 K=1,IPS
   40 IF (ISPTAB(I,N).EQ.IPT(K)) GOTO 30
      IPS=IPS+1
      IPT(IPS)=ISPTAB(I,N)
      KK=IPT(IPS)
      CALL SYMULT(SPS(1,1,IPS),TPS(1,IPS),S(1,1,KK),T(1,KK),
     1            SPS(1,1,1),TPS(1,1))
      I30 = 3
      CALL RECENT(TPS(1,IPS),I30)
C
C     IN ORDER TO GET THE POSITION CLOSEST TO THE ORIGIN
C     REDEFINE INDEPENDENT VARIABLES WHERE NECESSARY
C
      ATRAC=0.0
      DO 35 J=1,3
      ROW(J)=ABS(SPS(J,1,IPS))+ABS(SPS(J,2,IPS))+ABS(SPS(J,3,IPS))
      ATRAC=ATRAC+ABS(SPS(J,J,IPS))
   35 CONTINUE
      IF(INTEGER(ATRAC,0).EQ.(INTEGER(ATRAC1,0)-2)) THEN
        CALL TEHECU(SPS(1,1,1),SPS(1,1,IPS))
        ENDIF
      CALL REDEFI(SPS(1,1,IPS),TPS(1,IPS))
      IF (NLT.EQ.8) GOTO 30
      IF (GLAT.EQ.'H'.AND.SPS(3,3,IPS).LT.0.01) GOTO 30
      CALL SHCTOR(TPS(1,IPS),SPS(1,1,IPS),TU,NLT)
   30 CONTINUE
C
C     RESTORE THE COORDINATE TRIPPLE WITH...
C
      ISUM=0
C
C     14) NUL-MATRIX: SMALLEST TOTAL SHIFT
C
      TMSUM=0.0
      DO 45 I=1,3
      DO 45 II=1,3
   45 TMSUM=TMSUM+ABS(SPS(I,II,1))
      IF (TMSUM.GT.0.01) GOTO 60
      ISUM=IPS
      SUM=3.0
      DO 50 I=1,IPS
      ISS(I)=1
      SUM1=TPS(1,I)+TPS(2,I)+TPS(3,I)
   50 IF (SUM1.LT.(SUM-0.01)) SUM=SUM1
      DO 55 I=1,IPS
      SUM1=TPS(1,I)+TPS(2,I)+TPS(3,I)
      IF (INTEGER(SUM1,2).EQ.INTEGER(SUM,2)) GOTO 55
         ISUM=ISUM-1
         ISS(I)=-14
   55 CONTINUE
      IF (ISUM.EQ.1) GOTO 200
      GOTO 70
C
C     1) AT LEAST X POSITIVE WITH NO Y- OR Z-COMPONENT
C        AND FEWEST NUMBER OF NON-DIAGONAL ELEMENTS
C
   60 IOFF=3
      DO 80 I=1,IPS
      ISS(I)=1
      TEST=ABS(SPS(1,2,I))+ABS(SPS(1,3,I))
      IF (SPS(1,1,I).LT.-0.01.OR.TEST.GT.0.01) THEN
         ISS(I)=-1
      ELSE IF (SPS(1,1,I).GT.1.01) THEN
         ISS(I)=-1
      ELSE
         NOFF=INTEGER(ABS(SPS(2,1,I))+ABS(SPS(3,1,I))+ABS(SPS(3,2,I)),0)
         IF (NOFF.LT.IOFF) IOFF=NOFF
         ISUM=ISUM+1
         ENDIF
   80 CONTINUE
      DO 85 I=1,IPS
      IF (ISS(I).LT.0) GOTO 85
      NOFF=INTEGER(ABS(SPS(2,1,I))+ABS(SPS(3,1,I))+ABS(SPS(3,2,I)),0)
      IF (SPS(2,1,I).GT.1.01) NOFF=NOFF-1
      IF (NOFF.EQ.IOFF) GOTO 85
      ISUM=ISUM-1
      ISS(I)=-1
   85 CONTINUE
      IF (ISUM.EQ.1) GOTO 200
C
C     13) NOTHING LIKE X+1/2, Y+1/4, ...
C
   70 ISUM1=0
      DO 118 I=1,IPS
      IF (ISS(I).LT.0) GOTO 118
      DO 120 J=1,3
  120 IF (ABS(SPS(J,J,I)).GT.0.01.AND.TPS(J,I).GT.0.01) GOTO 122
      GOTO 118
  122 ISUM1=ISUM1+1
      ISS(I)=-13
  118 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.0) THEN
        DO 123 I=1,IPS
  123   IF(ISS(I).EQ.-13) ISS(I)=0
        ISUM=ISUM1
        ENDIF
      IF (ISUM.EQ.1) GOTO 200
C
C     2) Y .EQ. 0 OR 1
C
      ISUM1=0
      DO 90 I=1,IPS
      IF (ISS(I).LT.0) GOTO 90
      IF (SPS(2,2,I).GT.-0.01.AND.SPS(2,2,I).LT.1.01) GOTO 90
      ISUM1=ISUM1+1
      ISS(I)=-2
   90 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.0) THEN
         DO 95 I=1,IPS
   95    IF (ISS(I).EQ.-2) ISS(I)=0
         ISUM=ISUM1
         ENDIF
      IF (ISUM.EQ.1) GOTO 200
C
C     3) NO Z-COMPONENT IN Y
C
      ISUM1=0
      DO 100 I=1,IPS
      IF (ISS(I).LT.0) GOTO 100
      IF (SPS(2,3,I).EQ.0.0) GOTO 100
      ISUM1=ISUM1+1
      ISS(I)=-3
  100 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.0) THEN
         DO 105 I=1,IPS
  105    IF (ISS(I).EQ.-3) ISS(I)=0
         ISUM=ISUM1
         ENDIF
      IF (ISUM.EQ.1) GOTO 200
C
C     4) SMALLEST TRANSLATION IN X-DIRECTION
C
      TMIN=1.0
      ISUM1=0
      DO 110 I=1,IPS
      IF (ISS(I).LT.0) GOTO 110
      IF (TPS(1,I).LT.TMIN) TMIN=TPS(1,I)
  110 CONTINUE
      TMIN=TMIN+0.01
      DO 115 I=1,IPS
      IF (ISS(I).LT.0) GOTO 115
      IF (TPS(1,I).LT.TMIN) GOTO 115
      ISS(I)=-4
      ISUM1=ISUM1+1
  115 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.1) GOTO 200
C
C     5) Z .EQ. 0 OR 1
C
      ISUM1=0
      DO 135 I=1,IPS
      IF (ISS(I).LT.0) GOTO 135
      IF (SPS(3,3,I).GT.-0.01.AND.SPS(3,3,I).LT.1.01) GOTO 135
      ISUM1=ISUM1+1
      ISS(I)=-5
  135 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.0) THEN
         DO 140 I=1,IPS
  140    IF (ISS(I).EQ.-5) ISS(I)=0
         ISUM=ISUM1
         ENDIF
      IF (ISUM.EQ.1) GOTO 200
C
C     6) IF X (SPS(1,1,I)) .EQ. 0, NO X-COMPONENT OF Y
C
      ISUM1=0
      DO 138 I=1,IPS
      IF (ISS(I).LT.0) GOTO 138
      IF (SPS(1,1,I).NE.0.0) GOTO 138
      IF (SPS(2,1,I).EQ.0.0) GOTO 138
      ISS(I)=-6
      ISUM1=ISUM1+1
  138 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.0) THEN
         DO 142 I=1,IPS
  142    IF (ISS(I).EQ.-6) ISS(I)=0
         ISUM=ISUM1
         ENDIF
      IF (ISUM.EQ.1) GOTO 200
C
C     7) X-COMPONENT OF Y .LE. 1
C
      ISUM1=0
      DO 125 I=1,IPS
      IF (ISS(I).LT.0) GOTO 125
      IF (SPS(2,1,I).LT.1.01) GOTO 125
      ISUM1=ISUM1+1
      ISS(I)=-7
  125 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.0) THEN
         DO 130 I=1,IPS
  130    IF (ISS(I).EQ.-7) ISS(I)=0
         ISUM=ISUM1
         ENDIF
      IF (ISUM.EQ.1) GOTO 200
C
C     12) X-COMPONENT OF Y .GE. 0
C
      ISUM1=0
      DO 170 I=1,IPS
      IF (ISS(I).LT.0) GOTO 170
      IF (SPS(2,1,I).GT.-0.01) GOTO 170
      ISUM1=ISUM1+1
      ISS(I)=-12
  170 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.0) THEN
         DO 175 I=1,IPS
  175    IF(ISS(I).EQ.-12) ISS(I)=0
         ISUM=ISUM1
         ENDIF
      IF (ISUM.EQ.1) GOTO 200
C
C     8) SMALLEST TRANSLATION IN Y-DIRECTION
C
      TMIN=1.0
      DO 145 I=1,IPS
      IF (ISS(I).LT.0) GOTO 145
      IF (TPS(2,I).LT.TMIN) TMIN=TPS(2,I)
  145 CONTINUE
      TMIN=TMIN+0.01
      ISUM1=0
      DO 150 I=1,IPS
      IF (ISS(I).LT.0) GOTO 150
      IF (TPS(2,I).LT.TMIN) GOTO 150
      ISUM1=ISUM1+1
      ISS(I)=-8
  150 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.1) GOTO 200
C
C     9) SMALLEST TRANSLATION IN Z-DIRECTION
C
      TMIN=1.0
      DO 160 I=1,IPS
      IF (ISS(I).LT.0) GOTO 160
      IF (TPS(3,I).LT.TMIN) TMIN=TPS(3,I)
  160 CONTINUE
      TMIN=TMIN+0.01
      ISUM1=0
      DO 165 I=1,IPS
      IF (ISS(I).LT.0) GOTO 165
      IF (TPS(3,I).LT.TMIN) GOTO 165
      ISS(I)=-9
      ISUM1=ISUM1+1
  165 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.1) GOTO 200
C
C     10) IF X .EQ. 0, NO X-COMPONENT IN Z
C
      ISUM1=0
      DO 180 I=1,IPS
      IF (ISS(I).LT.0) GOTO 180
      IF (SPS(1,1,I).NE.0.0) GOTO 180
      IF (SPS(3,1,I).EQ.0.0) GOTO 180
      ISUM1=ISUM1+1
      ISS(I)=-10
  180 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.0) THEN
         DO 185 I=1,IPS
  185    IF (ISS(I).EQ.-10) ISS(I)=0
         ISUM=ISUM1
         ENDIF
      IF (ISUM.EQ.1) GOTO 200
C
C     11) NO NEGATIVE COMPONENTS IN Z
C
      ISUM1=0
      DO 187 I=1,IPS
      IF (ISS(I).LT.0) GOTO 187
      IF (SPS(3,1,I).GT.-0.01.AND.SPS(3,2,I).GT.-0.01) GOTO 187
      ISUM1=ISUM1+1
      ISS(I)=-11
  187 CONTINUE
      ISUM=ISUM-ISUM1
      IF (ISUM.EQ.0) THEN
         DO 189 I=1,IPS
  189    IF (ISS(I).EQ.-11) ISS(I)=0
         ISUM=ISUM1
         ENDIF
      IF (ISUM.EQ.1) GOTO 200
C
C     SEVERAL POSITIONS ARE EQUAL
C
      J=0
      DO 400 I=1,IPS
      IF (ISS(I).LT.0) GOTO 400
      IF (J.EQ.0) THEN
         J=I
         GOTO 400
         ENDIF
      IF (ICOMMA(SPS(1,1,J),SPS(1,1,I),3,2).EQ.0) THEN
         ISS(I)=-19
         ISUM=ISUM-1
         IF (ISUM.EQ.1) GOTO 200
         ENDIF
  400 CONTINUE
C
C     MAY BE: X,0,X  X,X,0  0,X,X  ETC.
C
      J=0
      DO 410 I=1,IPS
      IF (ISS(I).LT.0) GOTO 410
      IF (J.EQ.0) THEN
         J=I
         GOTO 410
         ENDIF
      DO 420 K=1,2
         DO 430 L=1,3
            DO 430 M=1,3
            N1=M+K
            IF(N1.GT.3) N1=N1-3
            SS(N1,L)=SPS(M,L,I)
  430       CONTINUE
         IF(ICOMMA(SPS(1,1,J),SS,3,2).EQ.0) THEN
            DO 440 L=1,3
               DO 440 M=1,3
               IF (SPS(L,M,J).LT.(SPS(L,M,I)-0.1)) THEN
                  M1=M+1
                  IF (M.EQ.3) M1=1
                  IF (SPS(M,M,J).GT.0.1.AND.SPS(M1,M,J).GT.0.1) THEN
                     ISS(I)=-18
                     ISUM=ISUM-1
                     IF (ISUM.EQ.1) GOTO 200
                     GOTO 420
                  ELSE IF (SPS(M,M,I).GT.0.1.AND.SPS(M1,M,I).GT.0.1)
     1               THEN
                     ISS(J)=-18
                     ISUM=ISUM-1
                     IF (ISUM.EQ.1) GOTO 200
                     J=I
                     GOTO 420
                  ELSE IF (SPS(M,M,J).LT.0.1.OR.SPS(M,M,I).LT.0.1) THEN
                     ISS(I)=-18
                     ISUM=ISUM-1
                     IF (ISUM.EQ.1) GOTO 200
                     ENDIF
                  ENDIF
  440          CONTINUE
            ENDIF
  420    CONTINUE
  410 CONTINUE
C
C     ERROR IN SORTING, NON-UNIQUE RESULT
C
      WRITE (NOUT,2000) ISUM
 2000 FORMAT (//,1X,'*****ERROR***SPESOR** NO UNIQUE SOLUTION,',I4,
     1   ' SEITZ MATRICIES LEFT',//,' *****CHOOSE YOUR OWN SOLUTION:')
      DO 190 I=1,IPS
      IF (ISS(I).LT.0) GOTO 190
      WRITE (NOUT,2010) I,(SPS(1,K,I),K=1,3),TPS(1,I),ISS(I),
     1  (SPS(2,K,I),K=1,3),TPS(2,I),IPT(I),(SPS(3,K,I),K=1,3),TPS(3,I)
 2010 FORMAT (/,1X,'NUMBER:',I4,10X,3F7.4,3X,F7.4,/,1X,'P-FLAG:',I4,
     1  10X,3F7.4,3X,F7.4,/,1X,'SYM-OP:',I4,10X,3F7.4,3X,F7.4)
      ISS(I)=-20
  190 CONTINUE
  195 WRITE (NOUT,2020)
 2020 FORMAT (/,1H$,'GIVE NUMBER OF MATRIX AS REPRESENTATIVE OF',
     1   ' WYKOFF GROUP: ')
      READ (NIN,1000,ERR=197) NSR
 1000 FORMAT (I5)
      IF (NSR.LT.0.OR.NSR.GT.IPS) THEN
         WRITE(NOUT,2040)
 2040    FORMAT(/,' ***** INPUT ERROR, TRY AGAIN *****',/)
         GOTO 195
         ENDIF
      ISS(NSR)=0
      GOTO 200
  197 WRITE(NOUT,2040)
      GOTO 195
C
C     WE FOUND THE GROUP REPRESENTATIVE
C
  200 IF (NDEBUG.NE.0) THEN
         DO 205 I=1,IPS
  205    WRITE (NDEBUG,2010) I,(SPS(1,K,I),K=1,3),TPS(1,I),ISS(I),
     1   (SPS(2,K,I),K=1,3),TPS(2,I),IPT(I),(SPS(3,K,I),K=1,3),TPS(3,I)
         STIME=SECOND(3)-STIME
         WRITE (NDEBUG,2030) STIME
 2030    FORMAT (/,' ++DEBUG++SPESOR++2030++ THIS SORT DONE IN ',
     1      F8.3,' SECONDS ')
         ENDIF
      DO 210 I=1,IPS
      IF (ISS(I).LT.0) GOTO 210
      SSIG=1.0
      TRACE=SPS(1,1,I)+SPS(2,2,I)+SPS(3,3,I)
      IF(TRACE.LT.-0.01) SSIG=-1.0
      DO 215 J=1,3
      TP(J,N)=TPS(J,I)
      DO 215 K=1,3
  215 SP(J,K,N)=SPS(J,K,I)*SSIG
      GOTO 10
  210 CONTINUE
   10 CONTINUE
C
C     SORT SPECIAL POSITIONS ACCORDING TO MULTIPLICITIES, AXIS
C     AND TRANSLATIONS
C
      STIME=SECOND(3)
      DO 250 I=2,NSP
      J=I
  260 J1=J-1
      IF (MULTSP(J).GE.MULTSP(J1)) THEN
         IF (MULTSP(J).GT.MULTSP(J1)) GOTO 300
         TRACE1=ABS(SP(1,1,J))+ABS(SP(2,2,J))+ABS(SP(3,3,J))
         TRACE2=ABS(SP(1,1,J1))+ABS(SP(2,2,J1))+ABS(SP(3,3,J1))
         IF (TRACE1.GT.(TRACE2+0.01)) GOTO 300
         IF (TRACE1.LT.(TRACE2-0.01)) GOTO 250
         OTR1=ABS(SP(2,1,J))+ABS(SP(3,1,J))+ABS(SP(3,2,J))
         OTR2=ABS(SP(2,1,J1))+ABS(SP(3,1,J1))+ABS(SP(3,2,J1))
         IF (OTR1.GT.(OTR2+0.01)) GOTO 300
         IF (OTR1.LT.(OTR2-0.01)) GOTO 250
         IF (TRACE1.GT.1.01) THEN
            IF (ABS(SP(3,3,J)).LT.(ABS(SP(3,3,J1))-0.01)) GOTO 300
            IF (ABS(SP(3,3,J)).GT.(ABS(SP(3,3,J1))+0.01)) GOTO 250
            IF (ABS(SP(2,2,J)).LT.(ABS(SP(2,2,J1))-0.01)) GOTO 300
            IF (ABS(SP(2,2,J)).GT.(ABS(SP(2,2,J1))+0.01)) GOTO 250
            IF (ABS(SP(1,1,J)).LT.(ABS(SP(1,1,J1))-0.01)) GOTO 300
            IF (ABS(SP(1,1,J)).GT.(ABS(SP(1,1,J1))+0.01)) GOTO 250
         ELSE
            IF (ABS(SP(3,3,J)).GT.(ABS(SP(3,3,J1))+0.01)) GOTO 300
            IF (ABS(SP(3,3,J)).LT.(ABS(SP(3,3,J1))-0.01)) GOTO 250
            IF (ABS(SP(2,2,J)).GT.(ABS(SP(2,2,J1))+0.01)) GOTO 300
            IF (ABS(SP(2,2,J)).LT.(ABS(SP(2,2,J1))-0.01)) GOTO 250
            IF (ABS(SP(1,1,J)).GT.(ABS(SP(1,1,J1))+0.01)) GOTO 300
            IF (ABS(SP(1,1,J)).LT.(ABS(SP(1,1,J1))-0.01)) GOTO 250
            ENDIF
         LR1=0
         LR2=0
         DO 290 M=1,3
            IF(TP(M,J).LT.0.01) LR1=LR1+1
  290       IF(TP(M,J1).LT.0.01) LR2=LR2+1
         IF(LR2.GT.LR1) GOTO 300
         SUM1=TP(1,J)+TP(2,J)+TP(3,J)
         SUM2=TP(1,J1)+TP(2,J1)+TP(3,J1)
         IF (SUM1.GT.(SUM2+0.01)) GOTO 300
         SUM1=SUM1-TP(3,J)
         SUM2=SUM2-TP(3,J1)
         IF (SUM1.GT.(SUM2+0.01)) GOTO 300
         IF (TP(1,J).GT.(TP(1,J1)+0.01)) GOTO 300
         ENDIF
      GOTO 250
  300 DO 310 K=1,NSYM
      ITEMP=ISPTAB(K,J1)
      ISPTAB(K,J1)=ISPTAB(K,J)
  310 ISPTAB(K,J)=ITEMP
      DO 330 K=1,3
      TEMP=TP(K,J1)
      TP(K,J1)=TP(K,J)
      TP(K,J)=TEMP
      DO 320 L=1,3
      TEMP=SP(K,L,J1)
      SP(K,L,J1)=SP(K,L,J)
  320 SP(K,L,J)=TEMP
  330 CONTINUE
      ITEMP=MULTSP(J1)
      MULTSP(J1)=MULTSP(J)
      MULTSP(J)=ITEMP
      DO 335 M=1,2
      GTEMP(M)=GSYSP(M,J1)
      GSYSP(M,J1)=GSYSP(M,J)
  335 GSYSP(M,J)=GTEMP(M)
      J=J1
      IF (J.LE.1) GOTO 250
      GOTO 260
  250 CONTINUE
C
C     PRINT STATEMENTS IN CASE NDEBUG .NE. 0
C
      IF (NDEBUG.EQ.0) RETURN
      STIME=SECOND(3)-STIME
      WRITE (NDEBUG,2200) STIME
 2200 FORMAT (//,' ++DEBUG++SPESOR++2200++ SPECIAL POSITIONS',
     1   ' SORTED IN ',F8.3,' SECONDS',//,1X,
     2   'MULTIPLICATION TABLE =ISPTAB= AT END OF =SPESOR= :',/)
      DO 360 I=1,NSP
  360 WRITE (NDEBUG,2210) (ISPTAB(J,I),J=1,NSYM)
 2210 FORMAT (2X,42I3,/,2X,6I3)
      WRITE (NDEBUG,2220) (MULTSP(I),I=1,NSP)
 2220 FORMAT (//,1X,'MULTIPLICITIES OF SPECIAL POSITIONS: ',//,
     1   1X,20I4)
      WRITE (NDEBUG,2230)
 2230 FORMAT (//,1X,'DEFINING SEITZ MATRICES OF SPECIAL POSITIONS:',
     1   //,1X,'NUMBER  ROTATIVE',13X,'TRANSLATIVE  PART',/)
      DO 340 I=1,NSP
  340 WRITE (NDEBUG,2240) I,((SP(K,L,I),L=1,3),TP(K,I),K=1,3)
 2240 FORMAT (1X,I4,3X,3F7.4,3X,F7.4,/,2(8X,3F7.4,3X,F7.4,/))
      WRITE (NDEBUG,2250) SECOND(3)
 2250 FORMAT (//,' ++DEBUG++SPESOR++2250++ END SPESOR AT',F8.3,/)
      RETURN
      END
C
C
C
C
      SUBROUTINE SPSYM
C     ----------------
C
C     WRITTEN BY D. ALTERMATT
C     MCMASTER UNIVERSITY, NOVEMBER 1984
C
C     CALLS INTEGER
C
C     THIS SUBROUTINE DERIVES THE POINT GROUP SYMMETRY FOR ALL
C     THE SPECIAL POSITONS STORED IN COMMON /SYMA/. A POINT GROUP
C     SYMBOL AS IN THE INTERNATIONAL TABLES FOR CRYSTALLOGRAPHY
C     IS CREATED AND STORED IN =GSYSP(2,I)=
C
C**** THE SPECIAL POSITION MULTIPLICATION TABLE =ISPTAB= MUST BE
C     CALCULATED BEFORE RUNNING =SPSYM= !!
C
      IMPLICIT CHARACTER*8 (G)
      IMPLICIT CHARACTER*80 (H)
C
C   SYSTEM COMMON BLOCKS
C
      COMMON /FILES/ NIN,NOUT,NDEBUG
      COMMON /SYM/   NSYM,MSYM,S(3,3,48)
      COMMON /GSYM/  GLAT,GCENT,HMSYM,HALL,GSYSP(2,30)
      COMMON /SYMA/  MULT(48,48),IT(3,48),NSP,MSP,ISPTAB(48,30)
C
C   =GGSYM= WILL CONTAIN THE SYMMETRIES IN
C   [100],[010],[001],[110],[-110],[111] AND [000]-DIRECTION
C
      DIMENSION GGSYM(7)
C
C   FIRST INTERPRET HALL SYMBOL TO FIND CRYSTAL SYSTEM
C   NCRYSY = 0 FOR TRIC. - ORTHO.
C            1 FOR TETRAGONAL
C            2 FOR HEX. PRIMITIVE
C            3 FOR HEX. R
C            4 FOR RHOMBOHEDRAL
C            5 FOR CUBIC
C
      NCRYSY=0
      DO 10 I=1,24
        K=ICHAR(HALL(I:I))
        IF(K.LT.ICHAR('A').OR.K.GT.ICHAR('R')) GOTO 10
        DO 20 J=I,24
          K=ICHAR(HALL(J:J))
          IF(K.LT.ICHAR('2').OR.K.GT.ICHAR('6')) GOTO 20
          IF(K.EQ.ICHAR('3')) THEN
            NCRYSY=2
            IF(HALL(J+1:J+1).EQ.'*') NCRYSY=4
            IF(GLAT.EQ.'H') NCRYSY=3
            GOTO 40
            ENDIF
          IF(K.EQ.ICHAR('4')) NCRYSY=1
          IF(K.EQ.ICHAR('6')) NCRYSY=2
          IF(GLAT.EQ.'H') NCRYSY=3
          IF(K.EQ.ICHAR('2').OR.K.EQ.ICHAR('4')) THEN
            DO 30 L=J,24
              K=ICHAR(HALL(L:L))
              IF(K.EQ.ICHAR('(')) GOTO 40
              IF(K.NE.ICHAR('3')) GOTO 30
              NCRYSY=5
              GOTO 40
   30         CONTINUE
            ENDIF
          GOTO 40
   20     CONTINUE
   10   CONTINUE
      IF(NDEBUG.NE.0) WRITE(NDEBUG,2030) NCRYSY
 2030 FORMAT(/,' ++DEBUG++SPSYM++2030++ NCRYSY =',I3)
C
C   FOR ALL SPECIAL POSITIONS DO...
C   WE ASSUME THAT THE SPECIAL POSITIONS ARE SORTED, STARTING WITH
C   THE GENERAL POSITION (I=1), ENDING WITH THE MOST SPECIAL POSI-
C   TION (I=NSP)
C
   40   DO 50 I=NSP,1,-1
        DO 60 J=1,7
   60   GGSYM(J)='        '
        NOSYM=1
C
C   SEARCH THROUGH =ISPTAB=
C
        DO 100 J=2,NSYM
          IF(ISPTAB(J,I).NE.1) GOTO 100
C
C   POSITION HAS THE SYMMETRY OF S(J)*S(1), GET IT
C
          NOSYM=0
          N=MULT(1,J)
          TRACE=S(1,1,N)+S(2,2,N)+S(3,3,N)
          ATRACE=ABS(S(1,1,N))+ABS(S(2,2,N))+ABS(S(3,3,N))
          OUTOT=S(1,2,N)+S(1,3,N)+S(2,1,N)+S(2,3,N)+S(3,1,N)+S(3,2,N)
          AOUTOT=ABS(S(1,2,N))+ABS(S(1,3,N))+ABS(S(2,1,N))+ABS(S(2,3,
     1           N))+ABS(S(3,1,N))+ABS(S(3,2,N))
          ITRAC=INTEGER(TRACE,0)
          IATRAC=INTEGER(ATRACE,0)
          IOUT=INTEGER(OUTOT,0)
          IAOUT=INTEGER(AOUTOT,0)
          IF(IABS(ITRAC).EQ.3) THEN
C   --      CENTRE OF SYMMETRY
            GGSYM(7)(1:2)='-1'
            GOTO 100
          ELSE IF(IATRAC.EQ.3.AND.IAOUT.EQ.0) THEN
            IF(ITRAC.EQ.-1) THEN
C   --        TWO-FOLD AXIS PARALLEL TO X, Y OR Z
              DO 110 K=1,3
  110         IF(S(K,K,N).GT.0.01) GGSYM(K)(4:4)='2'
              GOTO 100
            ELSE
C   --        MIRROR PLANE PERPENDICULAR TO X, Y OR Z
              DO 120 K=1,3
  120         IF(S(K,K,N).LT.-0.01) GGSYM(K)(8:8)='M'
              GOTO 100
              ENDIF
          ELSE IF(IATRAC.EQ.3) THEN
C   --      TWO-FOLD AXIS OR MIRROR PLANE IN [210], TREAT AS [110]
            IF(IOUT.GT.0) THEN
              DO 150 K=1,3
                IF(ITRAC.LT.-0.1) GGSYM(5)(4:4)='2'
  150           IF(ITRAC.GT.0.1)  GGSYM(4)(8:8)='M'
            ELSE
              DO 160 K=1,3
                IF(ITRAC.LT.-0.1) GGSYM(4)(4:4)='2'
  160           IF(ITRAC.GT.0.1)  GGSYM(5)(8:8)='M'
              ENDIF
            GOTO 100
          ELSE IF(IATRAC.EQ.2) THEN
            IF(ITRAC.EQ.0) THEN
C   --        3 OR 3BAR PARALLEL TO Z
              IF(GGSYM(3)(2:2).NE.'6') GGSYM(3)(2:2)='3'
              DO 170 K=1,3
                ROW=S(K,1,N)+S(K,2,N)+S(K,3,N)
                COL=S(1,K,N)+S(2,K,N)+S(3,K,N)
  170           IF(INTEGER(ROW,0).EQ.-1.AND.INTEGER(COL,0).EQ.-1)
     1            GGSYM(3)(1:1)='-'
              GOTO 100
            ELSE
C   --        6 OR 6BAR PARALLEL TO Z
              GGSYM(3)(2:2)='6'
              IF(ITRAC.LT.0) GGSYM(3)(1:1)='-'
              GOTO 100
              ENDIF
          ELSE IF(IATRAC.EQ.1) THEN
            IF(IABS(IOUT).EQ.IAOUT) THEN
C   --        TWO-FOLD AXIS OR MIRROR PLANE PARALLEL TO A DIAGONAL
              IF(IOUT.GT.0) THEN
                DO 130 K=1,3
                IF(S(K,K,N).LT.-0.1) GGSYM(4)(4:4)='2'
  130           IF(S(K,K,N).GT.0.1)  GGSYM(5)(8:8)='M'
                GOTO 100
              ELSE
                DO 135 K=1,3
                IF(S(K,K,N).LT.-0.1) GGSYM(5)(4:4)='2'
  135           IF(S(K,K,N).GT.0.1)  GGSYM(4)(8:8)='M'
                GOTO 100
                ENDIF
            ELSE
C   --        4 OR 4BAR PARALLEL TO Z
              DO 140 K=1,3
              K1=3
              IF(NCRYSY.EQ.5) K1=K
              IF(S(K,K,N).LT.-0.1) GGSYM(K1)(1:2)='-4'
  140         IF(S(K,K,N).GT.0.1)  GGSYM(K1)(2:2)='4'
              GOTO 100
              ENDIF
          ELSE
C   --      3 OR 3BAR PARALLEL TO BODY-DIAGONAL
            GGSYM(6)(2:2)='3'
            IF(IOUT.EQ.-3.OR.IOUT.EQ.1) GGSYM(6)(1:1)='-'
            ENDIF
  100     CONTINUE
C   --  IS IT THE GENERAL POSITION ?
        IF(NOSYM.EQ.1) GGSYM(7)(1:2)=' 1'
        IF(NDEBUG.NE.0) THEN
          WRITE(NDEBUG,2000) I,(ISPTAB(J,I),J=1,NSYM)
 2000     FORMAT(/,' ++DEBUG++SPSYM++2000++ SPEC.POS. NUMBER',I4,
     1      ' ISPTAB:',/,1X,48I2)
          WRITE(NDEBUG,2005) (GGSYM(K),K=1,7)
 2005     FORMAT(' GGSYM"  100 ""  010 ""  001 ""  110 ""',
     1    ' -110 ""  111 ""  000 "',/,6X,7A8)
          ENDIF
C
C   SO FAR, SO WELL. NOW GET RID OF EXESSIVE ELEMENTS
C
        IF(GGSYM(3)(1:4).EQ.'-6  ') GGSYM(3)(8:8)=' '
        ILEN=0
        DO 200 K=1,5
        IF(GGSYM(K)(8:8).EQ.'M') ILEN=ILEN+1
  200   IF(GGSYM(K)(1:2).NE.'  ') ILEN=ILEN+1
        IF(NCRYSY.NE.4.AND.ILEN.GE.2) THEN
          DO 210 K=1,5
  210     IF(GGSYM(K)(8:8).EQ.'M') GGSYM(K)(3:4)='  '
          ENDIF
C
C       UPDATE =NCRYSY= IN CASE HALL-SYMBOL SHOULD BE EMPTY
C
        IF(GGSYM(6).NE.'        ') THEN
          IF(NCRYSY.LT.4) NCRYSY=5
          IF(GGSYM(2)(2:2).EQ.'4') GGSYM(3)=GGSYM(2)
          GGSYM(2)='        '
          IF(NCRYSY.EQ.5) THEN
            IF(GGSYM(1)(2:2).EQ.'4') GGSYM(3)=GGSYM(1)
            GGSYM(1)='        '
            IF(GGSYM(4).NE.'        ') GGSYM(5)='        '
            ENDIF
          ENDIF
        DO 215 K=1,3
        IF(GGSYM(K)(1:2).NE.'  ') THEN
          GGSYM(K)(4:4)=' '
          IF(NCRYSY.LT.1.AND.GGSYM(K)(2:2).EQ.'4') NCRYSY=1
          IF(NCRYSY.LT.2.AND.GGSYM(K)(2:2).EQ.'3') NCRYSY=2
          IF(NCRYSY.LT.2.AND.GGSYM(K)(2:2).EQ.'6') NCRYSY=2
          ENDIF
        IF(GGSYM(K)(8:8).NE.' ') GGSYM(K)(1:1)=' '
  215   CONTINUE
        DO 220 K=1,6
  220   IF(GGSYM(K).NE.'        ') GGSYM(7)='        '
        DO 230 K=1,5
  230   IF(GGSYM(K)(1:4).NE.'    '.AND.GGSYM(K)(8:8).EQ.'M')
     1    GGSYM(K)(6:6)='/'
        IF(NCRYSY.EQ.5) THEN
          IF(GGSYM(1)(2:2).EQ.'4') THEN
            GDUM=GGSYM(3)
            GGSYM(3)=GGSYM(1)
            GGSYM(1)=GGSYM(2)
            GGSYM(2)=GDUM
            ENDIF
          IF(GGSYM(2)(2:2).EQ.'4') THEN
            GDUM=GGSYM(3)
            GGSYM(3)=GGSYM(2)
            GGSYM(2)=GGSYM(1)
            GGSYM(1)=GDUM
            ENDIF
          ENDIF
        IF(NDEBUG.NE.0) WRITE(NDEBUG,2005) (GGSYM(K),K=1,7)
C
C   ADD PERIODS, CONCATINATE, REMOVE BLANKS AND SAVE TO =GSYSP(2,I)=
C
        HLINE=' '
        IF(GGSYM(7).NE.' ') THEN
C   --    1 OR 1BAR ONLY
          HSYM=GGSYM(7)
          GOTO 310
          ENDIF
        IF(NCRYSY.EQ.0) THEN
C   --    SYMMETRY LESS THAN TETRAGONAL
          DO 300 K=1,3
  300     IF(GGSYM(K).EQ.' ') GGSYM(K)(8:8)='.'
          HSYM=GGSYM(1)//GGSYM(2)//GGSYM(3)//GGSYM(7)
          GOTO 310
          ENDIF
C   --    HIGHER SYMMETRIES
        DO 320 K=1,6
          IF(K.EQ.2) GOTO 320
          IF(GGSYM(K).EQ.' ') GGSYM(K)(8:8)='.'
  320     CONTINUE
          IF(NCRYSY.EQ.2.OR.NCRYSY.EQ.3) THEN
C   --      TRIGONAL OR HEXAGONAL
            IF(GGSYM(3)(1:2).NE.'  '.AND.GGSYM(4)(6:6).EQ.'/')
     1        GGSYM(4)(1:6)='      '
            IF(GGSYM(3)(1:2).NE.'  '.AND.GGSYM(5)(6:6).EQ.'/')
     1        GGSYM(5)(1:6)='      '
            IF(NCRYSY.EQ.3) THEN
              IF(GGSYM(4).NE.'        '.AND.GGSYM(5).EQ.'       .')
     1          GGSYM(5)='        '
              IF(GGSYM(5).NE.'        '.AND.GGSYM(4).EQ.'       .')
     1          GGSYM(4)='        '
              IF(GGSYM(4).EQ.GGSYM(5)) GGSYM(5)='        '
              ENDIF
            HSYM=GGSYM(3)//GGSYM(4)//GGSYM(5)//GGSYM(7)
            GOTO 310
            ENDIF
          IF(NCRYSY.EQ.1.OR.NCRYSY.GE.4) THEN
C   --      TETRAGONAL, RHOMBOHEDRAL OR CUBIC
            IF(NCRYSY.EQ.1) GGSYM(6)='        '
            IF(GGSYM(4).NE.'        '.AND.GGSYM(5).EQ.'       .')
     1        GGSYM(5)='        '
            IF(GGSYM(5).NE.'        '.AND.GGSYM(4).EQ.'       .')
     1        GGSYM(4)='        '
            IF(GGSYM(3)(1:2).NE.'  '.OR.NCRYSY.EQ.4) THEN
              GGSYM(2)='        '
              IF(GGSYM(4)(8:8).EQ.'M') GGSYM(5)='        '
              IF(GGSYM(5)(8:8).EQ.'M') GGSYM(4)='        '
              IF(GGSYM(4).EQ.GGSYM(5)) GGSYM(5)='        '
              ENDIF
            IF(GGSYM(1)(8:8).EQ.'.'.AND.GGSYM(2).NE.'        ')
     1        GGSYM(1)(8:8)=' '
            ENDIF
          IF(NCRYSY.EQ.4) THEN
C   --      RHOMBOHEDRAL
            IF(GGSYM(6)(1:4).NE.'    '.AND.GGSYM(4)(6:6).EQ.'/')
     1        GGSYM(4)(1:6)='      '
            IF(GGSYM(6)(1:4).NE.'    '.AND.GGSYM(5)(6:6).EQ.'/')
     1        GGSYM(5)(1:6)='      '
            HSYM=GGSYM(6)//GGSYM(4)//GGSYM(5)//GGSYM(7)
            GOTO 310
            ENDIF
          IF(NCRYSY.EQ.5) THEN
C   --      CUBIC
            IF(GGSYM(1)(8:8).EQ.'.') GGSYM(1)(8:8)=' '
            IF(GGSYM(3)(8:8).EQ.'M'.AND.GGSYM(6)(1:2).EQ.'-3')
     1        GGSYM(3)(1:6)='      '
            IF(GGSYM(4)(8:8).EQ.'M'.AND.GGSYM(6)(1:2).EQ.'-3')
     1        GGSYM(4)(1:6)='      '
            IF(GGSYM(5)(8:8).EQ.'M'.AND.GGSYM(6)(1:2).EQ.'-3')
     1        GGSYM(5)(1:6)='      '
            IF(GGSYM(4)(4:4).EQ.'2'.AND.GGSYM(5)(8:8).EQ.'M') THEN
              GGSYM(4)='       M'
              GGSYM(5)='   2    '
              ENDIF
            IF(GGSYM(1)(4:4).EQ.'2') THEN
              IF(GGSYM(2)(8:8).EQ.'M') THEN
                GGSYM(1)='       M'
                GGSYM(2)='   2    '
              ELSE IF(GGSYM(3)(8:8).EQ.'.') THEN
                GGSYM(3)='        '
              ELSE IF(GGSYM(3)(2:2).EQ.'4') THEN
                IF(GGSYM(4)(8:8).EQ.'.'.OR.GGSYM(5)(8:8).EQ.'.')
     1            GGSYM(1)='        '
                IF(GGSYM(3)(8:8).EQ.'M') THEN
                  GGSYM(3)(6:6)=' '
                  GGSYM(1)='        '
                  ENDIF
                ENDIF
              ENDIF
            ENDIF
        HSYM=GGSYM(3)//GGSYM(1)//GGSYM(2)//GGSYM(6)//GGSYM(4)//GGSYM(5)
     1    //GGSYM(7)
  310   IF(NDEBUG.NE.0) WRITE(NDEBUG,2010) HSYM
 2010   FORMAT(' ++2010++ HSYM: ',A80)
        J=1
        DO 250 K=1,48
          IF(HSYM(K:K).EQ.' ') GOTO 250
            HLINE(J:J)=HSYM(K:K)
            J=J+1
  250     CONTINUE
        GSYSP(2,I)=HLINE(1:8)
        IF(NDEBUG.NE.0) WRITE(NDEBUG,2020) HLINE,GSYSP(2,I)
 2020   FORMAT(' ++2020++HLINE: ',A80,/,' ++2020++ GSYSP(2,I): ',A8,/)
   50   CONTINUE
      RETURN
      END
C
C
C
C
      FUNCTION ITESTR(TA,TB,TU,NLT,SP)
C     ---------------
C
C   WRITTEN 84.5.11 BY I.D.BROWN
C   UPDATED 84.12.17 BY D.ALTERMATT
C
C   CALLS ICOMVE, REDEFI
C
C   EQUALS ZERO IF TA IS THE SAME AS ANY DERIVATIVE OF TB CREATED
C   WHEN THE SYMMETRY TRANSLATIONS OF TU ARE APPLIED TO ROWS OF TB.
C   ORIGIN OF TB IS REDEFINED WITH RESPECT TO MATRIX SP. TO PREVENT
C   THIS REDEFINITION SET ALL ELEMENTS OF SP TO 0.0
C
      DIMENSION TA(3),TB(3),TC(3),TU(3,NLT),SP(3,3)
      ITESTR=1
      DO 200 I=8,NLT
        DO 100 J=1,3
  100     TC(J)=TB(J)+TU(J,I)
        CALL REDEFI(SP,TC)
        IF(ICOMVE(TA,TC,3,2).EQ.0) THEN
          ITESTR=0
          RETURN
          ENDIF
  200   CONTINUE
      RETURN
      END
C
C
C
C
      SUBROUTINE RECENT(X,N)
C     ----------------------
C
C   WRITTEN 84.3.5 BY I D BROWN
C
C   BRINGS EACH ELEMENT OF X(N) INTO THE RANGE 0 TO 0.99 BY ADDING
C   OR SUBTRACTING 1
C
      COMMON/FILES/NIN,NOUT,NDEBUG
      DIMENSION X(N)
      DO 100 K=1,N
   10   IF(X(K).LE.1.0.AND.X(K).GE.0.0) GOTO 50
        IF(X(K).GE.1.) X(K)=X(K)-1.0
        IF(X(K).LT.0.) X(K)=X(K)+1.0
        GOTO 10
   50   IF(X(K).GT.0.9999) X(K)=0.0
  100   CONTINUE
      RETURN
      END
C
      SUBROUTINE LTRANS(NLT,T,TU)
C     -----------------
C
C   WRITTEN 84.2.8 BY I D BROWN
C
C   ADDS CENTERING TRANSLATIONS TO ARRAY TU
C
      DIMENSION T(3),TU(3,14)
      J=0
      NLT=NLT+1
      DO 50 I=1,3
        TU(I,NLT)=T(I)
        IF(T(I).LT.0.1) J=I
   50   CONTINUE
      IF(J.EQ.0) RETURN
C
C   IF ONE COMPONENT OF T IS ZERO ADD A UNIT TRANSLATION
C
      NLT=NLT+1
      DO 100 I=1,3
        TU(I,NLT)=T(I)
        IF(I.EQ.J) TU(I,NLT)=TU(I,NLT)+1.
  100   CONTINUE
      RETURN
      END
C
C
C
C
C
      SUBROUTINE XALTNO(XX,NPOS,NALT,NALTS,ICENT,NT,ZZ)
C
C     WRITTEN 85.7.26 BY I.D.BROWN
C
C     CALLS symop,recent,vminv,vplusv
C
C     This subroutine transforms xx to the altermatt position zz.
c     If NPOS.ne.0 only special position with multiplicity.eq.npos
c     will be searched, otherwise all will be searched and the correct
c     value of NPOS will be returned.
C     NALT is the Altermatt number found
C     NALTS is the symmetry operator to transform xx to zz
C     ICENT is the corresponding lattice centering operator
C     NT(3) is the subsequent lattice translation needed
c
c
*  *********** PROGRAMMERS NOTE *************************
*  THIS ROUTINE DOES NOT APPARENTLY RECOGNISE 0,Y,Z (E.G. 0,.09,.18)
*  AS X,Y,0 (POSTION 2) IN FM-3C
***********************************************************
C
      implicit character*8 (g)
      implicit character*80 (h)
c
c     system common blocks
c
      common /files/  nin,nout,ndebug,nread,nwrite,nspgp
      common /sym/    nsym,msym,s(3,3,48),t(3,48),nlt,tu(3,4)
      common /gsym/   glat,gcent,hmsym,hall,gsysp(2,30)
      common /syma/   mult(48,48),it(3,48),nsp,msp,isptab(48,30),
     1                  multsp(30),sp(3,3,30),tp(3,30)
       dimension xx(3),zz(3),nt(3),y(3),z(3),tx(3)
c
      nalt=0
c
c   check multiplicities and npos to screen possible special positons
c
      do 400 j=nsp,1,-1
        if(npos.eq.multsp(j).or.npos.eq.0) then
c         multiplicity of special position matches
          if(j.eq.1.and.npos.eq.0) then
c           the atom must be on a general position
              nalt=1
              nalts=1
              icent=1
              npos=multsp(1)
              do 100 k=1,3
                nt(k)=2
                zz(k)=xx(k)
  100           continue
              return
            endif
c
c   apply all symmetry operations and translations to xx to find a match
c   for sp/tp
c
          do 300 k=1,nsym
            call symop(y,s(1,1,k),t(1,k),xx)
            i30 = 3
            call recent(y,i30)
            do 250 k2=1,nlt
              call vminv(z,y,tu(1,k2))
              if(ialtst(sp(1,1,j),tp(1,j),z).eq.1)then
c             we have a match for the coordinates
                nalt=j
                if(npos.eq.0) npos=multsp(j)
                nalts=k
                icent=k2
                go to 500
                endif
  250         continue
  300       continue
          endif
  400   continue
      if(nalt.eq.0) return
c
c    Altermatt number (nalt), symmetry operator (nalts) and the lattice
c    centering translation (icent) have now been found.  We now need to
c    find the lattice translation and determine zz
c
  500 call symop(y,s(1,1,nalts),t(1,nalts),xx)
      call vplusv(z,y,tu(1,icent))
c
c   z should differ from the Altermatt coordinates only by a lattice
c   translation
c     First we generate the coordinates in the altermatt position
      call symop(zz,sp(1,1,nalt),tp(1,nalt),z)
      call recent(zz,3)
c     then find the lattice translation between zz and z
      call vminv(tx,zz,z)
      do 550 i=1,3
  550   nt(i)=integer(tx(i),0)+2
      return
      end
C
C
C
C
      SUBROUTINE NALTNO(XX,NPOS,NALT,NALTS,ICENT,NT)
C     ----------------------------------------------
C
C     WRITTEN BY I.D. BROWN, NOV. 22 1984
C
C     CALLS SYMULT, RECENT, IALTST, INTEGER, SYMOP, VMINV
C
C     THIS SUBROUTINE RETURNS IN =NALT= THE NUMBER OF THE SPECIAL POSITION
C     ('ALTERMATT-NUMBER') THAT MATCHES A GIVEN COORDINATE TRIPLET =XX=
C     OF MULTIPLICITY =NPOS=. IF =NPOS= IS 0, THEN ALL SPECIAL POSITIONS
C     ARE TESTED. =NALTS= WILL CONTAIN THE SYMMETRY OERATOR NUMBER FOR
C     THE TRANSFORMATION =SP(NALT)= TO =XX=, =ICENT= GIVES THE CENTRING
C     AND =NT(3)= THE SHIFTS TO NEIGHBOUR CELLS [(2,2,2) FOR NO SHIFT]
C
C     NALT=0 FOR NO MATCH
C
C**** RUN SPESOR BEFORE USING THIS SUBROUTINE.
C
      IMPLICIT CHARACTER*8 (G)
      IMPLICIT CHARACTER*80 (H)
C
C     SYSTEM COMMON BLOCKS
C
      COMMON /FILES/ NIN,NOUT,NDEBUG
      COMMON /SYM/   NSYM,MSYM,S(3,3,48),T(3,48),NLT,TA(3,4)
      COMMON /GSYM/  GLAT,GCENT,HMSYM,HALL,GSYSP(2,30)
      COMMON /SYMA/  MULT(48,48),IT(3,48),NSP,MSP,ISPTAB(48,30),
     1       MULTSP(30),SP(3,3,30),TP(3,30)
C
      DIMENSION SN(3,3),X(3),Y(3),XX(3),Z(3),TX(3),NT(3),TN(3),YY(3)
C
      NALT=0
C
C     LOOP THROUGH ALL SPECIAL POSITIONS, STARTING AT THE MOST
C     SPECIAL ONE
C
        DO 400 J=NSP,1,-1
          IF(NDEBUG.NE.0) WRITE(NDEBUG,2160) J,NPOS,MULTSP(J)
 2160     FORMAT(/,' ++DEBUG++NALTNO++2160++ J =',I3,
     1      ', NPOS(I) =',I3,', MULTSP(J) =',I3)
          IF(NPOS.EQ.MULTSP(J).OR.NPOS.EQ.0) THEN
C   --      MULTIPLICITIES OF ATOMIC AND SPECIAL POSITION MATCH
            IF(J.EQ.1.AND.NPOS.NE.0) THEN
C   --        GENERAL POSITION
              NALT=J
              NALTS=1
              ICENT=1
              DO 100 K=1,3
  100         NT(K)=2
              RETURN
              ENDIF
C
C     APPLY ALL SYMMETRY OPERATORS TO THE COORDINATES
C
            DO 300 K=1,NSYM
              CALL SYMULT(SN,TN,S(1,1,K),T(1,K),SP(1,1,J),TP(1,J))
              I30 = 3
              CALL RECENT(TN,I30)
              IF(NDEBUG.NE.0) WRITE(NDEBUG,2150) J,K,NPOS,MULTSP
     1          (J),((SN(J1,J2),J2=1,3),TN(J1),XX(J1),J1=1,3)
 2150         FORMAT(/' ++DEBUG++NALTNO++2150++ J =',I3,
     1          ',K =',I3,/,' NPOS(I) =',I3,', MULTSP(J) =',I3,
     2          ',      SN(3,3)     ,   TP(3),  XX(3)',3(/,30X,3F6.3
     3          ,3X,F6.3,3X,F6.3))
              DO 250 K2=1,NLT
                CALL VMINV(Z,XX,TA(1,K2))
                IF(IALTST(SN,TN,Z).EQ.1) THEN
                  ICENT=K2
                  NALT=J
                  NALTS=K
                  GOTO 500
                  ENDIF
  250           CONTINUE
  300         CONTINUE
            ENDIF
  400     CONTINUE
  500 IF(ICENT.NE.1) CALL VMINV(Z,XX,TA(1,ICENT))
      CALL SYMULT(SN,TN,S(1,1,NALTS),T(1,NALTS),SP(1,1,NALT),TP(1,NALT))
      CALL VMINV(Y,Z,TN)
      DO 530 I=1,3
        X(I)=0.0
        DO 520 J=1,3
          IF(ABS(SN(J,I)).GT.0.001) THEN
            X(I)=Y(J)/SN(J,I)
            GOTO 530
            ENDIF
  520     CONTINUE
  530   CONTINUE
      I30 = 3
      CALL RECENT(X,I30)
      CALL SYMOP(YY,SN,TN,X)
      CALL VMINV(TX,Z,YY)
      DO 550 I=1,3
  550   NT(I)=INTEGER(TX(I),0)+2
      RETURN
      END
C
C
C
C
      FUNCTION IALTST(SN,TM,X)
C
C     CALLS RECENT
C
C     SLAVE OF NALTNO
C
      DIMENSION SN(3,3),TM(3),X(3),Y(3),Z(3)
C
      IALTST=0
      DO 100 L=1,3
      Y(L)=0.0
  100 Z(L)=X(L)
      I30 = 3
      CALL RECENT(Z,I30)
      DO 200 L=1,3
        AROW=ABS(SN(L,1))+ABS(SN(L,2))+ABS(SN(L,3))
        IF(AROW.LT.0.001.AND.ABS(Z(L)-TM(L)).GT.0.001)
     1    GOTO 250
        IF(ABS(SN(L,1)).GT.0.001) Y(L)=(X(L)-TM(L))/SN(L,1)
  200   CONTINUE
      I30 = 3
      CALL RECENT(Y,I30)
      DO 300 L=1,3
        IF(ABS(SN(L,1)).LT.0.01) GOTO 300
        K1=L-1
        K2=L+1
        IF(K1.LT.1) K1=3
        IF(K2.GT.3) K2=1
        IF(ABS(SN(K1,1)).GT.0.001.AND.ABS(Y(K1)-Y(L)).GT.0.001)
     1    GOTO 250
        IF(ABS(SN(K2,1)).GT.0.001.AND.ABS(Y(K2)-Y(L)).GT.0.001)
     1    GOTO 250
  300   CONTINUE
      IALTST=1
  250 CONTINUE
      RETURN
      END