C          FILE GETSPEC
C               =======
C
C     V     V    AA     X     X
C     V     V   A  A     X   X
C      V   V   AAAAAA      X        -  VERSION
C       V V    A    A    X   X
C        V     A    A   X     X
C
C   THIS FILE CONTAINS THE FOLLOWING ROUTINES:
C
C  P GETSPEC        MAIN PROGRAM, INTERACTIVELY REQUESTS SPACE GROUP 
C                   SYMBOLS AND CALLS SUBROUTINES TO CALCULATE THE 
C                   SYMMETRY OPERATORS AND SPECIAL POSITIONS 
C  S SGHALL         CALCULATES SYMMETRY OPERATORS FROM HALL SYMBOL
C  S SYMD           DOUBLES NUMBER OF SYM.OPS. DUE TO CENTRE OF SYMMETRY
C  S SYMXYZ         PRINTS SYM.OPS. IN (X,Y,Z)-FORM
C  S GRMULT         CALCULATES THE POINT GROUP MULTIPLICATION TABLE
C  S SYMOP          DOES SYMMETRY TRANSFORMATIONS
C  S SYMULT         MULTIPLIES SYMMETRY OPERATORS
C  S WSCFS          CREATES AN OUTPUT FILE IN STANDARD CRYSTALLOGRAPHIC 
C                   FILE STRUCTURE (SCFS-84)
C  S INCREA         PRINTS OUT SYMMETRY TABLE
C  S CONXYZ         CONVERTS A SEITZ-MARIX TO (X,Y,Z)-FORM
C  F ICOMMA         COMPARES TWO (NXN)-MATRICES
C  F ICOMVE         COMPARES TWO (N)-VECTORS
C  F INTEGER        RETURNS INTEGER VALUE FOR REAL NUMBER, CORRECT ROUND-OFF
C  S MMINM          DIFFERENCE OF TWO MATRICES
C  S MXM            PRODUCT OF TWO MATRICES
C  S MXV            PRODUCT OF MATRIX * VECTOR
C  F SECOND         SYSTEM-DEPENDENT TIME CALLS
C  S VMINV          DIFFERENCE OF TWO VECTORS
C  S VPLUSV         SUM OF TWO VECTORS
C
C
c   GETSPEC calculates symmetry operators and special positions for any
c   given space group from the HALL space group symbol (which can be
c   obtained from the HERMANN-MAUGUIN symbol stored in the file SPGR.DAT).
c   The program is described by Altermatt and Brown in Acta Cryst. (1987)
c   A34, 125-130.
c
C   'GETSPEC' CALLS THE ROUTINES IN THE FILE 'GTSPLB'
C
C   'GETSPEC' OPENS THE EXISTING DATA FILE 'SPGR.DAT' IN DEFAULT DIRECTORY
C             CREATES FILE 'SCFS.DAT' AND, IF SPECIFIED, FILE 'DEBUG.DAT'
C             IN DEFAULT DIRECTORY
C
C   'SPGR.DAT' CONTAINS A CONCORDANCE BETWEEN HERMANN-MAUGUIN AND HALL
C             SPACE GROUP SYMBOLS
C
C
      PROGRAM GETSPEC
C
C   ***** INTERACTIVE VAX-VERSION *****
C
C   THIS PROGRAM WILL PROMPT FOR A HERMANN-MAUGUIN SYMBOL AS USED IN THE 
C   INORGANIC CRYSTAL STRUCTURE DATABASE, THE HALL SPACE GROUP SYMBOL AND
C   THE SPACE GROUP NUMBER. ONLY ONE OF THE TWO SYMBOLS IS NEEDED AND THE
C   SPACE GROUP NUMBER IS OPTIONAL. THE PROGRAM WILL THEN SEARCH THROUGH 
C   THE DATAFILE 'SPGR.DAT' TO GET THE EQUIVALENT OTHER SYMBOL AND SPACE 
C   GROUP NUMBER (IF NOT ALREADY GIVEN), CALCULATE THE SYMMETRY OPERATORS 
C   AND SPECIAL POSITIONS, SORT, PRINT AND STORE (ON FILE 'SCFS.DAT') ALL 
C   THE DATA.
C
C   THE FOLLOWING ROUTINES ARE CALLED:
C
C       GETSPEC:  SECOND, SGHALL, GRMULT, WSCFS, CONXYZ
C       SPELIB:   SPOST, SPESOR, SPMULT, SPSYM
C
C
C   THE FINAL ORDERING OF THE SPECIAL POSITIONS IS THAT DESCRIBED BY THE
C   AUTHOR IN ACTA CRYST (1985)
C
C
C    NOTES ON THE NAMES OF ARRAYS USED IN MCMASTER 
C         CRYSTALLOGRAHIC PROGRAMS
C
C    NOTE:  Matrices are given as M(row,column)
C           default integer names are used
C           implicit character*8 (g)
C           implicit character*80 (h)
C    FILES
C    -----
C    nin         input file number
C    nout        output file number
C    ndebug      debug output file number(set to 0)
C    nread       input data file number
C    nwrite      output data file number
C    nspgr       input space group file
C
C    SYM
C    ---
C    nsym        number of symmetry operators stored
C    msym        dimension of symmetry arrays (=48)
C    s(3,3,48)   rotation matrix
C    t(3,48)     translation matrix
C    nt          multiplicity of lattice
C    tu(3,4)     translations due to lattice centrings (including (0,0,0))
C
C    GSYM
C    ----
C    glat        lattice type
C    gcent       A/C acentric/centric
C    hmsym       Hermann-Mauguin space group symbol
C    hall        Hall space group symbol
C    gsysp(2,30) Wyckoff symbol, site symmetry
C
C    SYMA
C    ----
C    mult(48,48) symmetry multiplication table
C    it(3,48)    symmetry translations in integer form (*12)
C    nsp         number of special positions stored
C    msp         dimensions of special position arrays (=30)
C    isptab(48,30) equivalent symmetry operators for a given special
C                  position
C    multsp(30)  multiplicity of special postions
C    sp(3,3,30)  rotation matirix for special position operator
C    tp(3,30)    translation vector for secial position operator
C
C    REFERENCE:  ACTA CRYST.
C
C    WRITTEN IN JUNE 1985 BY DANIEL ALTERMATT, MCMASTER UNIVERSITY, HAMILTON
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)
C
C     START TIMER
C
      TIME = SECOND(1)
C
C     GREET USER AND SET DEFAULT VALUES IN COMMON BLOCKS
C
      NIN = 5
      NOUT = 6
      WRITE(NOUT,2000)
 2000 FORMAT('1',////,10x,' Welcome to program GETSPEC !',//,
     1  10x,' We will prompt you for a Hermann-Mauguin and a Hall',/,
     2  10x,' Space Group Symbol and then derive all symmetry',/,
     3  10x,' operators and special positions for the given space',/,
     4  10x,' group. ',//,
     5  10x,' You may enter either one or both symbols. If given,',/,
     6  10x,' the Hall-Symbol will be used for calculations.',///)
C-----------------
C     IF YOU WOULD LIKE A DEBUG OUTPUT, YOU MAY SET 'NDEBUG = 8' AND 
C     REMOVE THE 'C' FROM THE FOLLOWING LINE. THIS WILL WRITE DEBUG 
C     STATEMENTS TO 'DEBUG.DAT'.
C     ***** FOR CUBIC GROUPS THE SIZE OF 'DEBUG.DAT' MAY GET ENORMOUS.
C
      NDEBUG = 0
C     OPEN (8,FILE='DEBUG.DAT',STATUS='NEW')
C-----------------
      NREAD = NIN
      NWRITE = 9
      OPEN (9,FILE='SCFS.DAT',STATUS='NEW',CARRIAGECONTROL='LIST')
      NSPGR = 10
      OPEN (10,FILE='SPGR.DAT',STATUS='UNKNOWN')
C
C     STARTING POINT FOR RE-RUNS, GET CURRENT RUNTIME
C
   10 TIME = SECOND(3)
C
      NSYM = 0
      MSYM = 48
      NT = 1
      TU(1,1) = 0.0
      TU(2,1) = 0.0
      TU(3,1) = 0.0
C
      GLAT = ' '
      GCENT = ' '
      HMSYM = ' '
      HALL = ' '
C
      NSP = 0
      MSP = 30
C
C     CHECK IF 'SPGR.DAT' EXISTS IN CURRENT DIRECTORY
C
      REWIND(NSPGR)
      READ(NSPGR,1000,END=900) HLINE
 1000 FORMAT(A80)
      REWIND(NSPGR)
   15 WRITE(NOUT,2010)
 2010 FORMAT(/,'$ Enter Herman-Mauguin Space Group Symbol (or <RTN>): ')
      READ(NIN,1000,ERR=20,END=20) HMSYM
      IF (HMSYM.EQ.'<RTN>'.OR.HMSYM.EQ.'<rtn>') GOTO 20
      GOTO 30
   20 HMSYM = ' '
   30 WRITE(NOUT,2020)
      HMS = HMSYM
 2020 FORMAT(/,'$ Enter Hall Space Group Symbol (or <RTN>): ')
      READ(NIN,1000,ERR=40,END=40) HALL
      IF (HALL.EQ.'<RTN>'.OR.HALL.EQ.'<rtn>') GOTO 40
      GOTO 50
   40 HALL = ' '
   50 IF (HALL.EQ.' '.AND.HMSYM.EQ.' ') THEN
         WRITE(NOUT,2030)
 2030    FORMAT(//'  Come on Buddy, enter at least one legal symbol',
     1      ' !!!',///)
         GOTO 15
         ENDIF
C
C     CONVERT INPUT TO CAPITAL LETTERS
C
      IF (HALL.NE.' ') THEN
         HMS = ' '
         DO 45 I=1,15
            IF (ICHAR(HALL(I:I)).LT.ICHAR('a').OR.ICHAR(HALL(I:I))
     1         .GT.ICHAR('z')) GOTO 45
            INP = ICHAR(HALL(I:I))-ICHAR('a')+ICHAR('A')
            HALL(I:I) = CHAR(INP)
   45       CONTINUE
      ELSE
         DO 55 I=1,24
            IF (ICHAR(HMSYM(I:I)).LT.ICHAR('a').OR.ICHAR(HMSYM(I:I))
     1         .GT.ICHAR('z')) GOTO 55
            INP = ICHAR(HMSYM(I:I))-ICHAR('a')+ICHAR('A')
            HMSYM(I:I) = CHAR(INP)
   55       CONTINUE
         ENDIF
      WRITE(NOUT,2040)
 2040 FORMAT(/,'  Enter Space Group Number (as in International',
     1   ' Tables',/,'$ for Crystallography) or 0 (zero): ')
      READ(NIN,1010,ERR=60,END=60) NGROUP
 1010 FORMAT(I5)
      IF (NGROUP.LT.1.OR.NGROUP.GT.230) GOTO 60
      GOTO 70
   60 NGROUP = 0
C
C     MAKE SURE THAT THE SYMBOLS ARE LEFT JUSTIFIED
C
   70 IF (HMS .NE. ' ') THEN
         DO 75 I=1,24
   75       IF (HMSYM(I:I).NE.' ') GOTO 80
   80    II = 24-I+1
         IF (II.EQ.24) GOTO 100
         HMSYM(1:II) = HMSYM(I:24)
         HMSYM(II+1:24) = ' '
      ELSE
         DO 85 I=1,15
   85       IF (HALL(I:I).NE.' ') GOTO 90
   90    II = 15-I+1
         IF (II.EQ.15) GOTO 100
         HALL(1:II) = HALL(I:15)
         HALL(II+1:15) = ' '
         ENDIF
C
C     READ 'HALL' AND/OR 'NGROUP' FROM 'NSPGR'
C
  100 READ(NSPGR,1000,END=910) HLINE
      IF (HALL.EQ.' ') THEN
         IF (HMSYM(1:24).NE.HLINE(6:29)) GOTO 100
         HALL(1:15) = HLINE(30:44)
      ELSE
         IF (HALL(1:15).NE.HLINE(30:44)) GOTO 100
         HMSYM(1:24) = HLINE(6:29)
         ENDIF
      READ(HLINE(45:49),'(BZ,I5)') NGRP
      IF (NGRP.EQ.0) GOTO 110
      IF (NGROUP.NE.0.AND.NGROUP.NE.NGRP) THEN
         WRITE(NOUT,2050) NGROUP,NGRP
 2050    FORMAT(//,' ***** ERROR: Space Group Number given is:',i5,
     1   //,13x,' Space Group Number found is:',i5,
     2   //,13x,' Second number taken !!',//)
         ENDIF
      NGROUP = NGRP
C
C     CALCULATE THE SYMMETRY OPERATORS
C
  110 WRITE(NOUT,2060)
 2060 FORMAT(//,'  ----- Calculating the Symmetry Operators ...')
      CALL SGHALL
      IER = 0
      CALL GRMULT(IER)
      IF (IER.NE.0) GOTO 930
C
C     CALCULATE AND SORT THE SPECIAL POSITIONS
C
      WRITE(NOUT,2070)
 2070 FORMAT(//,'  ----- Calculating the Special Positions ...')
      IER = 0
      CALL SPOST(IER,NDEBUG)
      IF (IER.NE.0) GOTO 920
      WRITE(NOUT,2080)
 2080 FORMAT(//,'  ----- Sorting the Special Positions ...')
      CALL SPESOR
      CALL SPMULT
      CALL SPSYM
C
C     OK, ALL DONE. PRINT AND STORE RESULTS
C
      CALL WSCFS
      WRITE(NOUT,2100) NGROUP,HMSYM(1:15),HALL(1:15)
 2100 FORMAT(20(/),'  Space Group Number             ',i5,/,
     1             '  Hermann-Mauguin Symbol         ',a15,/,
     2             '  Hall Symbol                    ',a15,///,
     3  '$ Hit <RETURN> to continue: ')
      READ(NIN,1020,ERR=120,END=120) GDUM
 1020 FORMAT(A1)
  120 WRITE(NOUT,2110)
 2110 FORMAT(///,'  Symmetry Operators:',/)
      LINC = 0
      DO 130 I=1,NSYM
         LINC = LINC+1
         CALL CONXYZ(S(1,1,I),T(1,I),HBUF)
         WRITE(NOUT,2120) I,HBUF
 2120    FORMAT(1X,I5,5X,A30)
         IF (LINC.EQ.20.OR.I.EQ.NSYM) THEN
            WRITE(NOUT,2130)
 2130       FORMAT(//,'$ Hit <RETURN> to continue: ')
            READ(NIN,1020,ERR=140,END=140) GDUM
  140       IF (I.NE.NSYM) WRITE(NOUT,2110)
            LINC = 0
            ENDIF
  130    CONTINUE
      WRITE(NOUT,2140)
 2140 FORMAT(///,'  Listing of Special Positions:',//,
     1   '  Nr.   Coordinates          Multiplicity    Site Symmetry',/)
      LINC = 0
      DO 150 I=1,NSP
         LINC = LINC+1
         CALL CONXYZ(SP(1,1,I),TP(1,I),HBUF)
         WRITE(NOUT,2150) I,HBUF,MULTSP(I),GSYSP(2,I)
 2150    FORMAT(1X,I4,3X,A21,5X,I3,11X,A8)
         IF (LINC.EQ.18.OR.I.EQ.NSP) THEN
            WRITE(NOUT,2130)
            READ(NIN,1020,ERR=160,END=160) GDUM
  160       IF (I.NE.NSP) WRITE(NOUT,2140)
            LINC = 0
            ENDIF
  150    CONTINUE
      TIME=SECOND(3)-TIME
      WRITE(NOUT,2155) TIME
 2155 FORMAT(//,'  Elapsed CPU-Time: ',F7.2,' Seconds')
  170 WRITE(NOUT,2160)
 2160 FORMAT(/////,'$ Another Space Group to be done (Y/N): ')
      READ(NIN,1020,ERR=180,END=180) GYN
      IF (GYN.EQ.'Y'.OR.GYN.EQ.'y') GOTO 10
      IF (GYN.EQ.'N'.OR.GYN.EQ.'n') GOTO 999
  180 WRITE(NOUT,2170)
 2170 FORMAT(//,'  ****** INPUT ERROR ******* TRY AGAIN ******',//,
     1          '  only YES (yes,Y,y) or NO (no,N,n) are valid')
      GOTO 170
  900 WRITE(NOUT,2200) NSPGR
 2200 FORMAT(///,'  ****** ERROR: File',i3,' empty ******',//,
     1   '  ''SPGR.DAT'' must be in your current Directory !!',////)
      GOTO 999
  910 IF (HALL.NE.' ') THEN
         IF(NGROUP.LE.1.OR.NGROUP.GT.230) NGROUP = 0
         GOTO 110
         ENDIF
      WRITE(NOUT,2210) NSPGR
 2210 FORMAT(///,'  ****** ERROR: End-of-File in unit',i3,' ******')
      IF (HMSYM.NE.' ') WRITE(NOUT,2220) 'Hermann-Mauguin',HMSYM
      IF (HALL.NE.' ')  WRITE(NOUT,2220) '           Hall',HALL
 2220 FORMAT(//,2X,A15,' Symbol ',A15,' not found in ''SPGR.DAT'' ',
     1   ////)
      GOTO 170
  920 WRITE(NOUT,2230) IER
 2230 FORMAT(///,'  ****** ERROR: Malfunction of Routine SPOST ******',
     1   /,'         Error number',i3,//,'         Now you''re really',
     2   ' screwed !!',//,'         Try to DEBUG by yourself or call',
     3   ' the author...',////)
      GOTO 999
  930 WRITE(NOUT,2240)
 2240 FORMAT(///,'  ****** ERROR: Your Space Group Symbol is probably',
     1   ' incomplete',////)
      GOTO 170
  999 CLOSE(NSPGR)
      CLOSE(NWRITE)
      WRITE(NOUT,2300)
 2300 FORMAT(/,'  Your current Directory now contains a file SCFS.DAT',
     1   /,'  where above information is stored in Standard Crystal-',
     2   /,'  lographic File Structure (SCFS).',////)
      STOP
      END
C
C
C
C
      SUBROUTINE SGHALL
*
*
* WRITTEN BY DAVID R. MOSSCROP (7901822)
* MAY-JUNE ; 1982.
*
* FROM AN ARTICLE AND PROGRAM: NEW SPACE GROUP NAMES.
*      BY HALL 1982  REF:ACTA CRYST .A36, 517-523.
*
* CALLS SYMD, SYMXYZ
*
* THIS PROCEDURE USES THE HALL SPACE GROUP NAME TO GENERATE 
* THE CORPSONDING MATRIX TRANSFORMATIONS.
C
C   MODIFIED 84.2.14 BY I D BROWN TO ALLOW ORIGIN SHIFT
C   UPDATED BY D. ALTERMATT
*
      IMPLICIT CHARACTER*8 (G)
      IMPLICIT CHARACTER*80 (H)
C
      INTEGER OLD
      INTEGER ROT(12,20)
      INTEGER SMAT(16,50)
      INTEGER SVEC(3)
      DIMENSION ISHIFT(3),GCHR(30)
      INTEGER TAB(3,9)
      INTEGER VEC(3,8)
      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
*
* VARIABLE DICTIONARY: 
*
* GCENT IS CENTROSYMMETRIC FLAG (N OR C)
* GLAT IS LATTICE FLAG (P,R,A,B,C,I,F). 
* IAXS IS CURRENT AXIS POINTER TO TAB.
* ISGN IS MATRIX SIGN FLAG.
* IROT IS CURRENT ORDER POINTER TO TAB. 
* JAXS IS PREVIOUS AXIS POINTER TO TAB. 
* JROT IS PREVIOUS ORDER POINTER TO TAB.
* NIN IS INPUT DEVICE NUMBER. 
* NOUT IS OUTPUT DEVICE NUMBER.
* NSYM IS NUMBER OF EQUIVALENT POSITIONS.
* NGEN IS NUMBER OF GENERATED SEITZ MATRICES.
* NROT IS CURRENT ROTATION ORDER (1,2,3,4,6).
* NEW IS POSITION OF MATRIX PRODUCT.
* OLD IS MATRIX COUNT BEFORE INNER LOOP.
* ROT(12,20) IS TABLE OF ROTATION MATRICES.
* SMAT(16,50) IS SEITZ MATRIX ARRAY.
* SVEC(3) IS TEMPORARY T-VECTOR ARRAY.
* ISHIFT IS THE ORIGIN SHIFT
* TAB(3,9) POINTS TO ROT(AXIS,ORDER).
* VEC(3,8) IS THE T-VETORS FOR A,B,C,N,U,V,W,D.
*
******************************************************************
* LOADING INITIAL DATA
*
      DATA ROT/
     1 1,0,0,0, 0,1,0,0, 0,0,1,0,
     2 1,0,0,0, 0,-1,0,0, 0,0,-1,0,
     3 -1,0,0,0, 0,1,0,0, 0,0,-1,0,
     4 -1,0,0,0, 0,-1,0,0, 0,0,1,0,
     5 1,0,0,0, 0,0,1,0, 0,-1,-1,0,
     6 -1,0,-1,0, 0,1,0,0, 1,0,0,0,
     7 0,1,0,0, -1,-1,0,0, 0,0,1,0,
     8 1,0,0,0, 0,0,1,0, 0,-1,0,0,
     9 0,0,-1,0, 0,1,0,0, 1,0,0,0,
     1 0,1,0,0, -1,0,0,0, 0,0,1,0,
     2 1,0,0,0, 0,1,1,0, 0,-1,0,0,
     3 0,0,-1,0, 0,1,0,0, 1,0,1,0,
     4 1,1,0,0, -1,0,0,0, 0,0,1,0,
     5 0,0,1,0, 1,0,0,0, 0,1,0,0,
     6 -1,0,0,0, 0,0,-1,0, 0,-1,0,0,
     7 0,0,-1,0, 0,-1,0,0, -1,0,0,0,
     8 0,-1,0,0, -1,0,0,0, 0,0,-1,0,
     9 -1,0,0,0, 0,0,1,0, 0,1,0,0,
     1 0,0,1,0, 0,-1,0,0, 1,0,0,0,  0,1,0,0, 1,0,0,0, 0,0,-1,0/
*
      DATA GCHR/'1','2','3','4','5','6','P','R','A','B','C','I','F',
     1'N','U','V','W','D',' ','-','X','Y','Z','*','''','"','(','0','H',
     2'H'/
*
      DATA TAB/1,1,1,2,3,4,5,6,7,8,9,10,0,0,0,11,12,13,                 
     1  14,14,14,15,16,17,18,19,20/
*
      DATA VEC/6,0,0, 0,6,0, 0,0,6, 6,6,6,                              
     1  3,0,0, 0,3,0, 0,0,3, 3,3,3/
*
* VARIABLE INITIALIZATION.
*
      JCENT=0
      GCENT=GCHR(14) 
      GLAT=' '
      NGEN=1
      NSYM=1
      ISGN=12
      NLT=1
*
      NROT=0
      LLTCODE=0
      IROT=0
      JROT=0
      IAXS=0
      JAXS=0
*
******************************************************************
* THE PROGRAM: SGHALL
*
* ARRAY INITIALIZATION.
*
*
* ZERO VECTOR ACCUMULATOR.
*
      DO 1 I=1,3
      SVEC(I)=0
      TU(I,NLT)=0.0
    1 CONTINUE
*
* SET THE (4,4) ELEMENT TO 1.0.
*
      DO 2 I=1,50
      SMAT(16,I)=12 
      SMAT(4,I)=0
      SMAT(8,I)=0
      SMAT(12,I)=0
    2 CONTINUE
*
* SET ORIGIN SHIFT TO ZERO
*
      DO 30 I=1,3
   30   ISHIFT(I)=0 
*
* SET IDENTITY ROTATIONS.
*
      DO 3 I=1,12
      SMAT(I,1)=ROT(I,1)*12
    3 CONTINUE
*
* SET IDENTITY TRANSLATIONS.
*
      DO 4 I=13,16
      SMAT(I,1)=ROT(I-10,1)*12
    4 CONTINUE
*
* CHECKING ELEMENT BY ELEMENT AND SET VARIABLES AS NEEDED.
*
      DO 66 N=1,15
      DO 6 M=1,29
      IF (HALL(N:N).EQ.GCHR(M)) GO TO 606
    6 CONTINUE
*
* ERROR CONDITION CHARACTER NOT IDENTIFIED.
*
      WRITE(NOUT,1670) HALL(N:N),N
 1670 FORMAT(' ++ERROR++SGHALL++1670++ ILLEGAL CHARACTER: ',A1,
     1       ' IN COLUMN ',I1)
*
* SKIP INACTIVE BLANKS.
*
  606 IF ((M.EQ.19).AND.(NROT.EQ.0)) GO TO 66
C
C  BLANK INITIATES CONSTRUCTION OF SYMMETRY GENERATOR
C
      IF ( M.EQ.19 ) GO TO 1011
C
C PARENTHESIS INITIATES READING OF ORIGIN SHIFT
C
      IF(M.EQ.27) GO TO 40
*
*
* SET GCENT AND GLAT, DEPENDENT ON FIRST ELEMENT(S).
*
      IF ( LLTCODE .NE. 0) GO TO 1000
*
* FIRST NEGATIVE SIGN INDICATES CENTRE. 
*
      IF (M .NE. 20) GO TO 1001
      JCENT=1
      GCENT=GCHR(11)
      GO TO 66
*
* ALLOWED CHARACTERS: P,R,A,B,C,I,F TO SET LTCODE.
*
 1001 IF ((M .GE. 7).AND.(M .LE.13)) GO TO 2000
      WRITE(NOUT,1600) GCHR(M)
 1600 FORMAT(' ++ERROR++SGHALL++1600++ UNRECOGNISED LATTICE TYPE ',
     1 A1)
      GO TO 150
 2000 CONTINUE
      LLTCODE=M-6
      GLAT=GCHR(M)
*
* SPACE GROUPS WITH NAMES 'R' TYPE ARE ASSIGNED
* GLAT 'H' FOR A HEXAGONAL SETTING. THIS IS A
* LOCAL CONVENTION. 
*
      IF ( M .EQ. 8 ) GLAT=GCHR(30)
*
* ADD CENTERING TRANSLATIONS TO =TU=
*
      IF(GLAT.EQ.'P') GOTO 660
      IF(GLAT.EQ.'I') THEN
        NLT=NLT+1
        DO 60 I=1,3
   60   TU(I,NLT)=0.5
        GOTO 660
        ENDIF
      IF(GLAT.EQ.'A'.OR.GLAT.EQ.'F') THEN
        NLT=NLT+1
        TU(1,NLT)=0.0
        TU(2,NLT)=0.5
        TU(3,NLT)=0.5
        IF(GLAT.EQ.'A') GOTO 660
        ENDIF
      IF(GLAT.EQ.'B'.OR.GLAT.EQ.'F') THEN
        NLT=NLT+1
        TU(1,NLT)=0.5
        TU(2,NLT)=0.0
        TU(3,NLT)=0.5
        IF(GLAT.EQ.'B') GOTO 660
        ENDIF
      IF(GLAT.EQ.'C'.OR.GLAT.EQ.'F') THEN
        NLT=NLT+1
        TU(1,NLT)=0.5
        TU(2,NLT)=0.5
        TU(3,NLT)=0.0
        GOTO 660
        ENDIF
      IF(GLAT.EQ.'H') THEN
        NLT=NLT+1
        TU(1,NLT)=0.3333333
        TU(2,NLT)=0.6666667
        TU(3,NLT)=0.6666667
        NLT=NLT+1
        TU(1,NLT)=0.6666667
        TU(2,NLT)=0.3333333
        TU(3,NLT)=0.3333333
        ENDIF
  660 IF(NDEBUG.NE.0) WRITE(NDEBUG,2020) ((TU(I1,I2),I1=1,3),I2=1,NLT)
 2020 FORMAT(/,' ++DEBUG++SGHALL++2020++  TU(3,NLT)',4(/,26X,3F5.2))
      GOTO 66
*
* IF CHARACTER IS A NEGATIVE SIGN, SET ISGN TO NEGATIVE,
* OTHERWISE POSITIVE. THE LATER NEGATIVE SIGNS INDICATE
* INVERSIONS.
*
 1000 IF (M.NE.20) GO TO 999
      ISGN=-12
      GO TO 66
*
* ROTATION MUST FOLLOW LATTICE TYPE.
*
*
* IF ROTATION AXIS ALREADY SET GO TO 1003.
*
  999 IF (NROT .NE. 0) GO TO 1003
      IF (M .LE.6) GO TO 2011 
      WRITE(NOUT,2010) GCHR(M)
 2010 FORMAT(' ++ERROR++SGHALL++2010++ ILLEGAL ROTATION CHARACTER '
     1 ,A1)
      GO TO 150
 2011 CONTINUE
*
* SET ROTATION AXIS IN NROT.
*
      IF ((M .EQ. 1) .AND. (ISGN .LT.0)) NSYM=NSYM*2
      NROT=M
      IROT=M
      IF ((M.NE.1).OR.(ISGN.LT.0)) NGEN=NGEN+1
      IF(NGEN.LE.MSYM) GO TO 1005
      WRITE(NOUT,1501)
 1501 FORMAT(' ++ERROR++SGHALL++1501++ NGEN GREATER THAN MSYM.')
      NSYM=0
      RETURN
*
*********** SELECT AXIS EXPLICILY **************
*
* IF CHARACTER DOES NOT SET AXIS GO TO 1005.
*
 1003 IF (M .LE. 20) GO TO 1009
      IAXS=M-20
      GO TO 1004
*
* SET AXIS IMPLICITLY.
*
*
*  IF DEFAULT AXIS NOT REQUIRED GO TO 1004.
*
 1005 CONTINUE
*
*  FIRST AXIS ASSUMED Z.
*
      IF ((NGEN .NE.2).AND.(IROT .NE.1)) GO TO 1006
      IAXS=3
      GO TO 1004
*
* SECOND AXIS ASSUMED X IF NROT=2 AND PREVIOUS AXIS WAS 2 OR 4.
*
 1006 IF ((NROT .NE.2) .OR. (JROT .NE.2)) GO TO 1007
      IAXS=1
      GO TO 1004
*
* AXIS AFTER 3,6 ASSUMED AS X-Y.
*
 1007 IF (( NROT .NE.2) .OR.(JROT .NE.3))  GO TO 1008
      IAXS=5
      GO TO 1004
 1008 IF (NROT .EQ.3) GO TO 2500
      WRITE(NOUT,1639)
 1639 FORMAT(' ++ERROR++SGHALL++1639++ UNABLE TO ASSIGN IMPLICIT ',
     1 'AXIS.')
      GO TO 150
 2500 IAXS=4
*
*
* FOR NON-PRINCIPAL AXES.
*
 1004 IF (IAXS .LE.3) GO TO 66
      IROT=IAXS+3
      IAXS=JAXS
      IF((IAXS.EQ.0).AND.(IROT.EQ.7))IAXS=1
      GO TO 66
*
********** SPECIFY TRANSLATION ***********
*
*
* ONLY FOR X,Y,Z DIRECTIONS.
*
 1009 IF(M .GT. 6) GO TO 1010 
      IF (IROT .GT.M) GO TO 2510
      WRITE(NOUT,1660) GCHR(M)
 1660 FORMAT(' ++ERROR++SGHALL++1660++ ILLEGAL SCREW AXIS SYMBOL ',
     1 A1)
      GO TO 150
 2510 CONTINUE
*
* CALCULATE TRANSLATIONS IN 12 THS.
*
      SVEC(IAXS)=SVEC(IAXS)+M*12/NROT
*
* TEST FOR A,B,C,N,U,V,W,D.
*
 1010 IF(((M.LE.8).OR.(M.GE.12)).AND.((M.LE.13).OR.(M.GE.19)))GOTO1011
      IF (M .GE. 12) GO TO 1111
      I=M-8
      GO TO 1112
 1111 I=M-10
 1112 DO 7 J=1,3
      SVEC(J)=SVEC(J)+VEC(J,I)
    7 CONTINUE
      GO TO 66
*
* CONSTRUCT GENERATOR WHEN BLANK IS REACHED.
*
 1011 IF (M .NE.19) GO TO 66
*
* POINT TO ROTATION MATIX.
*
      K=TAB(IAXS,IROT)
      DO 8 I=1,12
      SMAT(I,NGEN)=ROT(I,K)*ISGN
    8 CONTINUE
      DO 9 I=13,15
*
* DEPOSIT TRANSLATON BASE 12. 
*
      SMAT(I,NGEN)=MOD(SVEC(I-12)+120,12)
    9 CONTINUE
*
* CALCULATE THE NUMBER OF EQUIVALENT POSITIONS.
*
      NSYM=NSYM*NROT
      IF (NROT .GT. 3) NROT=NROT/2
*
* RESET ORDER FLAGS AND AXIS FLAGS.
*
      JROT=NROT
      NROT=0
      IROT=0
      JAXS=IAXS
      IAXS=0
      ISGN=12
*
* ZERO VECTOR ACCUMULATOR.
*
      DO 10 I=1,3
      SVEC(I)=0
   10 CONTINUE
*
      IF(HALL(N:N+1).EQ.'  ') GOTO 666
      GO TO 66
C
C   READ ORIGIN SHIFT
C
   40 IAX=0
      NN=N
   41 ISIGN=1
      IAX=IAX+1
      IF(IAX.GT.3) GO TO 666
   42 NN=NN+1
      DO 45 M=1,29
   45 IF(HALL(NN:NN).EQ.GCHR(M)) GO TO 47
      WRITE(NOUT,1670) HALL(NN:NN),NN
      GO TO 150
   47 IF(M.EQ.20) THEN
        ISIGN=-1
        GO TO 42
        ENDIF
      IF(M.EQ.28) THEN
        ISHIFT(IAX)=0
        GO TO 41
        ENDIF
      IF(M.LE.6) THEN
        ISHIFT(IAX)=ISIGN*M
        GO TO 41
        ENDIF
      IF(M.EQ.29) GO TO 666
      WRITE(NOUT,200) HALL(NN:NN),NN
      NSYM=0
      RETURN
   66 CONTINUE
  666 CONTINUE
      IF(NDEBUG.NE.0) THEN
        WRITE(NDEBUG,1666) GLAT,NROT,ISHIFT
        DO 667 I=1,NGEN
          WRITE(NDEBUG,1667) (SMAT(J,I),J=1,16)
  667     CONTINUE
 1666   FORMAT(' ++DEBUG++SGHALL++1666++ LATCODE=',A1,' NROT=',I3,
     1   ' ISHIFT=',3I4)
 1667   FORMAT(10X,16I4)
        ENDIF
*
* END OF INPUT DATA (HALL SPACE GROUP NAME) EVALUATION LOOPS.
*
* CHECK FOR ERROR CONDITION.
*
      IF ((NROT .LE.0).AND.(LLTCODE .NE.0)) GO TO 151
  150 CONTINUE
*
* ERROR CONDITION.
*
      WRITE(NOUT,200) HALL(N:N),N,NROT,LLTCODE,GLAT,GCENT
  200 FORMAT (' ++DEBUG++SGHALL++200++  ILLEGAL CONSTRUCTION AT',
     1       ' CHARACTER ',A1,' IN COL.',I3,/,25X,'NROT =',I3,
     1       ', LLTCODE =',I3,', GLAT = ',A1,', GCENT = ',A1,/)
      NSYM=0
      RETURN
  151 CONTINUE
*
*
*GENERATE MATRICES BY PROGRESSIVE BINARY PRODUCTS.
*
*
*
* SKIP GENERATOR IF P1 OR -P1.
*
      IF (NGEN .LE.1) GO TO 117
      M=2 
  723 IF(M.LE.NGEN) GO TO 273 
      GO TO 105
*
* SAVE CURRENT MATRIX COUNT.
*
  273 OLD=NGEN
      N=M 
  724 IF(N.LE.NGEN) GO TO 274 
      GO TO 11
*
* POINT TO GENERATED MATRIX IN SMAT.
*
  274 NEW=NGEN+1
*
      DO 13 I=1,3
      DO 131 J=1,13,4
  131 SMAT(I+J-1,NEW)=((SMAT(I,M)*SMAT(J,N))+(SMAT(I+4,M)*SMAT(J+1,N))+ 
     1   (SMAT(I+8,M)*SMAT(J+2,N))+(SMAT(I+12,M)*SMAT(J+3,N)))/12
      SMAT(I+12,NEW)=MOD(SMAT(I+12,NEW)+120,12)
   13 CONTINUE
*
* IS MATRIX UNIQUE? 
*
      LFLAG=1
      DO 14 L=1,NGEN
      DO 103 I=1,12 
  103 IF (SMAT(I,L).NE.SMAT(I,NEW)) GO TO 1674
      IF (LLTCODE .NE. 1) GO TO 251
      DO 15 II=13,16
      IF (SMAT(II,L) .EQ. SMAT(II,NEW)) GO TO 15
*
* ERROR CONDITION FOUND.
*
      WRITE (NOUT,300)
  300 FORMAT('++SGHALL++ERROR: MATRICES RELATED BY PURE TRANSLATION.')
      NSYM=0
      RETURN
   15 CONTINUE
*
      GO TO 251
*
 1674 DO 16 J=1,12
   16 IF (SMAT(J,L) .NE. SMAT(J,NEW)) GO TO 20
*
* ERROR CONDITION NOT AVOIDED ( IE NOT SKIPPED OVER).
*
      WRITE(NOUT,400)
  400 FORMAT('++SGHALL++ERROR: SOME MATRICES RELATED BY INVERSION.')
      NSYM=0
      RETURN
*
   20 CONTINUE
      LFLAG=LFLAG+1 
   14 CONTINUE
  251 IF (LFLAG .GT. NGEN) NGEN=NEW
      N=N+1
      GO TO 724
   11 IF ((NGEN.EQ.OLD).AND.(NGEN.GE. NSYM)) GO TO 105
      M=M+1
      GO TO 723
  105 CONTINUE
*
*
* SORT ON ZERO PATTERN.
*
      DO 17 N=2,NGEN
      DO 18 L=N,NGEN
      DO 181 K=1,12 
*
* EXIST ON ZERO MISMATCH.
*
  181 IF((SMAT(K,N-1).EQ.0).AND.(SMAT(K,L).NE.0))GO TO 18
      GO TO 1171
   18 CONTINUE
      GO TO 17
*
*
* EXCHANGE SORTED ELEMENTS.
*
 1171 DO 19 K=1,16
      I=SMAT(K,L)
      SMAT(K,L)=SMAT(K,N)
      SMAT(K,N)=I
   19 CONTINUE
*
   17 CONTINUE
  117 CONTINUE
*
* ERROR CHECK OF NGEN RELATIVE TO NSYM. 
*
      IF (NGEN .LE. NSYM) GO TO 1013
      WRITE(NOUT,500) NGEN,NSYM
  500  FORMAT('++SGHALL++ERROR: ',I5,' MATRICES GENERATED,',I5,
     1       ' EXPECTED')
      NSYM=0
*
* THE FOLLOWING LOOPS ENTERS THE VALUES OF SMAT(16,48),
* WHICH IS THE ARRAY USED IN THIS ROUTINE TO STORE THE
* ROTATION AND TRANSLATION MATRICES, INTO TWO COMPARABLE
* ARRAYS S(3,3,48) AND T(3,48) WHICH ARE USED IN OTHER
* ROUTIES EXTERNAL TO SGHALL. 
*
 1013 DO 7775 M=1,48
      DO 7775 J=1,3 
      DO 7775 I=1,3 
      K=4*(J-1)+I
*
* REAL DIVISION BY ZERO IS NEEDED TO CONVERT THE 12 THS
* USED IN SGHALL INTO ROTATIONS AND TRANSLATIONS
* UNDERSTANDABLE TO ROUTINE WHICH USE S(3,3,48) AND T(3,48).
* THE VALUES IN S AND T ARE NOW 1 OR DECIMAL FRACTIONS.
*
      S(I,J,M)=SMAT(K,M)/12.0 
      T(I,M)=SMAT(I+12,M)/12.0
 7775 CONTINUE
C
C  ADD ORIGIN SHIFT 
C
      DO 7780 M=1,48
        DO 7777 I=1,3
          IF(S(I,I,M).LT.-0.001) THEN
            T(I,M)=T(I,M)+ISHIFT(I)/6.0 
            IF(T(I,M).GT.0.99) T(I,M)=T(I,M)-1.0
            IF(T(I,M).LT.-0.01) T(I,M)=T(I,M)+1.0 
            ENDIF
 7777     CONTINUE
 7780   CONTINUE
      IF(GCENT.EQ.GCHR(11)) CALL SYMD(S,T,NSYM,MSYM)
      IF(NDEBUG.NE.0) THEN
        NOUT1=NOUT
        NOUT=NDEBUG
        CALL SYMXYZ (NSYM,S,T)
        NOUT=NOUT1
        ENDIF
      RETURN
      END 
C
C
C
C
      SUBROUTINE SYMD(S,T,NSYM,MSYM)
C
C  *****SUBROUTINE*****           SYMD
C                                 ----
C
C WRITTEN 80.6.1 BY DOUG THIERS
C
C SUBROUTINE TO DOUBLE THE NUMBER OF SYM.OP.S FOR CENTROSYMETRIC S.G.S
C
      DIMENSION S(3,3,MSYM),T(3,MSYM)
      COMMON/FILES/NIN,NOUT,NDEBUG,NREAD,NWRITE,NSPGP
      DO 5 I=1,NSYM 
      IF(NSYM+I.LE.MSYM) GOTO 10
      WRITE(NOUT,100)
  100 FORMAT(40H ERROR IN SUBROUTINE:SYMD.MORE SYM OP.S
     1   ,40HCALCULATED THAN ALLOWED(48).                             )
      RETURN
   10 DO 15 J=1,3
      DO 25 K=1,3
   25 S(J,K,NSYM+I)=-S(J,K,I) 
   15 T(J,NSYM+I)=T(J,I)
    5 CONTINUE
      NSYM=2*NSYM
      RETURN
C**  THIS PROGRAM VALID ON FTN4 AND FTN5 **
      END 
C
C
C
C
      SUBROUTINE SYMXYZ(NSYM,S,T)
*
* WRITTEN BY DAVID R. MOSSCROP (7901822)
* MAY-JUNE ; 1982.
C 84.6.12 D. ALTERMATT
C
* FROM AN ARTICLE AND PROGRAM: NEW SPACE GROUP NAMES.
*      BY HALL 1982  REF:ACTA CRYST .A36, 517-523.
*
* CALLS INTEGER
*
* THIS SUBROUTINE PRINTS OUT THE SPACE GROUP TRANSFORMATIONS
* IN FAMILIAR X,Y,Z COORDINATES.
*
*
* TRANS IS A ONE DIMENSIONAL ARRAY OF 30 ELEMENTS 
* IN TRANS ONE TRANSFORMATION IS STORED BEFORE
* BEING PRINTED OUT IN XYZ FORM.
*
      IMPLICIT CHARACTER*8 (G)
      IMPLICIT CHARACTER*80 (H)
*
      INTEGER TRANS(30),INUM(7),CHR(26) 
      INTEGER SMAT(16,48)
      INTEGER BLNK,COMA,OVER,PLUS,SUBT
      DIMENSION S(3,3,NSYM), T(3,NSYM)
      COMMON /FILES/ NIN,NOUT 
      COMMON /GSYM/ GLAT,GCENT,HMSYM,HALL
      DATA BLNK,COMA,PLUS,SUBT,OVER/1H ,1H,,1H+,1H-,1H//
      DATA INUM/1H1,1H2,1H3,1H5,1H3,1H4,1H6/
      DATA CHR/1H1,1H2,1H3,1H4,1H5,1H6,1HP,1HR,1HA,1HB,1HC,1HI,1HF,     
     1  1HN,1HU,1HV,1HW,1HD,1H ,1H-,1HX,1HY,1HZ,1H*,1H',1H"/
*
* IN THIS ROUTINE NGEN IS EQUAL TO NSYM.
*
      IF (NSYM .GT. 0) GO TO 1021
      WRITE(NOUT,*)' NOTHING TO PRINT.' 
      RETURN
*
* EXTERNALLY THE ROTATIONS AND TRANSLATIONS ARE STORED AND
* MANIPULATED IN ARRAYS S(3,3,48) AND T(3,48) RESPECTIVELY. 
* IN SYMXYZ ( THIS ROUTINE ) THE ARRAY SMAT(16,48) SERVES
* AS STORAGE FOR BOTH THE ROTATIONS AND TRANSLATIONS.
*
 1021 DO 7775 M=1,NSYM
      DO 7775 J=1,3 
      DO 7775 I=1,3 
      K=4*(J-1)+I
      SMAT(K,M)=INTEGER((S(I,J,M)*12),0)
*
* MULTIPLICATION BY 12 IS NECESSARY TO CHANGE THE DECIMAL
* FRACTIONS OF S(3,3,48) AND T(3,48) INTO THE 12 THS
* WHICH THIS ROUTINE REQUIRES.
* 12 THS ARE USED BECAUSE THE TRANSLATION FRACTIONS
* (1/4, 1/3, 3/4, 1/2, 1/6, 5/6 ETC) ALL ARE NICE, IE EVEN
* MULTIPLIES OF 12. 
*
      SMAT(I+12,M)=INTEGER((T(I,M)*12),0)
 7775 CONTINUE
*
* PRINT HALL SPACE GROUP SYMBOL.
*
      WRITE(NOUT,100)   HALL
  100 FORMAT(' HALL SPACE GROUP SYMBOL = ',A15)
*
      WRITE(NOUT,105)
  105 FORMAT('    ')
      WRITE(NOUT,110)
  110 FORMAT(' SYMMETRY OPERATORS ')
      DO 20 M=1,NSYM
*
* PRESET ARRAY TRANS TO BLANKS.
*
      DO 21 I=1,30
      TRANS(I)=BLNK 
   21 CONTINUE
      N=1 
*
* EXTRACT TRANSLATIONAL ELEMENTS.
*
      DO 22 L=1,3
      J=SMAT(L+12,M)
C
C SPECIAL POSITION WITH COORDINATE(S) '0'
C
      NSUM=IABS(SMAT(4*L-3,M))+IABS(SMAT(4*L-2,M))+IABS(SMAT(4*L-1,M))
     1     +IABS(J) 
      IF(NSUM.EQ.0) THEN
         TRANS(N)=1H0
         N=N+1
         ENDIF
      IF (J.LE.0) GO TO 1017
*
* SET THE NUMERATOR AS ONE (1).
*
      IF (J.GE.8) GO TO 1014
      TRANS(N)=INUM(1)
      GO TO 1015
*
* SET NUMERATOR AS 2, 3 OR 5. 
*
 1014 TRANS(N)=INUM(J-6)
*
* INSERT SLASH INTO FRACTION. 
*
 1015 TRANS(N+1)=OVER
*
* POINT TO DENOMINATOR CODE.
*
      I=IABS(6-J)
      IF (J.NE.6) GO TO 1016
*
* SET DENOMINATOR AS TWO (2)
*
      TRANS(N+2)=INUM(2)
      GO TO 1020
*
* SET DENOMINTATOR AS 3, 4 OR 6
 1016 TRANS(N+2)=INUM(I+3)
 1020 N=N+3
 1017 DO 23 K=1,3
      J=L+(K-1)*4
      J=SMAT(J,M)
      IF (J.EQ.0) GO TO 23
      IF (J.LE.0) GO TO 1018
*
* INSERT PLUS SIGN. 
*
      TRANS(N)=PLUS 
      GO TO 1019
*
* INSERT NEGATIVE SIGN.
*
 1018 TRANS(N)=SUBT 
*
* INSERT X, Y, OR Z.
*
 1019 TRANS(N+1)=CHR(K+20)
      N=N+2
   23 CONTINUE
*
* SEPARATE WITH COMMA.
*
      TRANS(N)=COMA 
      N=N+1
   22 CONTINUE
*
* REMOVE THE LAST COMMA FROM THE OUTPUT.
*
      TRANS(N-1)=BLNK
*
* LIST SYMMETRY AS X,Y,Z.
*
      WRITE(NOUT,700)M,TRANS
  700 FORMAT(3X,I3,5X,30A1)
   20 CONTINUE
      RETURN
      END 
C
C
C
C
      SUBROUTINE GRMULT(IER)
C     -----------------
C
C   WRITTEN 84.1.13 BY I.D.BROWN
C
C   CALLS MXM AND ICOMMA
C
C***** CALLS TIME-FUNCTION =SECOND(I)=
C
C   THIS ROUTINE CALCULATES THE POINT GROUP MULTIPLICATION TABLE "MULT"
C   FROM THE ARRAY "S" IN /SYM/
C
      COMMON /FILES/ NIN,NOUT,NDEBUG
      COMMON /SYM/   NSYM,MSYM,S(3,3,48)
      COMMON /SYMA/  MULT(48,48),IT(3,48),NSP,MSP,ISPTAB(48,30)
C
      DIMENSION SS(3,3)
C
      IER=0
      DO 200 I=1,NSYM
        DO 150 J=1,NSYM
          CALL MXM(SS,S(1,1,J),S(1,1,I),3)
          DO 100 K=1,NSYM
            IF(ICOMMA(S(1,1,K),SS,3,2).EQ.0) THEN 
              MULT(I,J)=K
              GO TO 150
              ENDIF 
  100       CONTINUE
          WRITE(NOUT,1000) I,J
 1000     FORMAT(' ++ERROR++GRMULT++1000++ NO MATCH FOR PRODUCT OF'
     +     ,' SYMMETRY OPERATORS ',2I3) 
          IF(NDEBUG.NE.0) WRITE(NDEBUG,1001) SS
 1001     FORMAT(' ++DEBUG++GRMULT++1001++ ',9F6.2)
          IER=1
          RETURN
  150    CONTINUE
  200 CONTINUE
      IF(NDEBUG.NE.0) THEN
        WRITE(NDEBUG,1200) SECOND(3)
 1200   FORMAT(' ++DEBUG++GRMULT++1200++(',F6.2,') GROUP ', 
     + 'MULTIPLICATION TABLE')
        DO 220 I=1,NSYM
  220     WRITE(NDEBUG,1201) (MULT(I,J),J=1,NSYM) 
 1201     FORMAT(5X,40I3,/,10X,8I3)
        ENDIF
      RETURN
      END 
C
C
C
C
      SUBROUTINE SYMOP(Y,S,T,X)
C     ----------------
C
C   WRITTEN 84.2.8 BY I D BROWN
C
C   APPLIES THE SEITZ MATRIX S:T TO X TO GIVE Y
C
      DIMENSION X(3),Y(3),S(3,3),T(3)
      DO 50 I=1,3
        Y(I)=0.
        DO 20 J=1,3 
  20      Y(I)=Y(I)+S(I,J)*X(J)
        Y(I)=Y(I)+T(I)
   50   CONTINUE
      RETURN
      END 
C
C
C
C
      SUBROUTINE SYMULT(S3,T3,S2,T2,S1,T1)
C     -----------------
C
C   WRITTEN 84.2.8 BY I D BROWN
C
C   CALLS MXM, MXV, VPLUSV
C
C   CALCULATES (S/T)3=(S/T)2*(S/T)1
C
      DIMENSION S1(3,3),S2(3,3),S3(3,3),T1(3),T2(3),T3(3),T4(3)
      CALL MXM(S3,S2,S1,3)
      CALL MXV(T4,S2,T1)
      CALL VPLUSV(T3,T4,T2)
      RETURN
      END 
*
*
*
*
      SUBROUTINE WSCFS
*
* JULY 20-24 82 BY DAVID R. MOSSCROP / 84.6.5 D. ALTERMATT
* / 86.7.11 N. HEINIG /
*
* REDUCED TO SYMMETRY INFORMATIONS FOR =GETSPEC=
*
* CALLS INTEGER, INCREA
*
* THIS SUBROUTINE PRINTS THE COMMON BLOCK CONTENTS IN STANDARD
* FORMAT, AS DESCRIBED IN " THE STANDARD CRYSTALLOGRAPHIC FILE
* STRUCTURE" BY I. D. BROWN.(ACTA CRYST (1983) A39, 216-224).
*
*
      IMPLICIT CHARACTER*8 (G)
      IMPLICIT CHARACTER*80 (H)
*
      COMMON /FILES/ NIN,NOUT,NDEBUG,NREAD,NWRITE,NSPGP
      COMMON /SYM/   NSYM,MSYM,S(3,3,48),T(3,48)
      COMMON /GSYM/  GLAT,GCENT,HMSYM,HALL,GSYSP(2,30)
      COMMON /SYMA/  MULT(48,48),IT(3,48),NSPE,MSPE,ISPTAB(48,30),
     1               MUSP(30),SP(3,3,30),TSP(3,30)
*
* TITLE SECTION
*
      WRITE(NWRITE,1000) HALL,HALL
 1000 FORMAT(//,'TITLE   ',/,'* Symmetry Operators and Special ',
     1  'Positions for ',A15,5X,A8,/)
*
* CHECK FOR SPACE GROUP SYMBOL(S), EITHER HALL AND/OR HERMANN-MAUGUIN.
* FIRST CHECK IF BOTH ARE PRESENT.
*
      IF(HALL(1:7).NE.'       '.OR. HMSYM(1:7).NE.'       '
     *   .OR. GLAT.NE.' ' .OR. GCENT.NE.' ') THEN 
        WRITE(NWRITE, 1040)
 1040   FORMAT('SG NAME')
*
* AT LEAST ONE PRESENT
*
        IF(HALL(1:8).NE.'        ') WRITE(NWRITE,1060) HALL 
 1060   FORMAT(1X,'HALL',5X,A24)
        IF(HMSYM(1:8).NE.'        ') WRITE(NWRITE,1065) HMSYM
 1065   FORMAT(1X,'HERM',5X,A24)
        IF(GLAT.NE.' ' .OR. GCENT.NE.' ') WRITE(NWRITE,1070)GCENT,GLAT
 1070   FORMAT(1X,'LATT',5X,A1,A1)
        WRITE(NWRITE,1050)
 1050   FORMAT('*EOS',/)
	ENDIF
*
* CHECK FOR SYMMETRY OPERATORS.
*
      IF(NSYM.NE.0) THEN 
        WRITE(NWRITE,1080)
 1080   FORMAT('SYMMETRY')
        CALL INCREA
*
* CHECK FOR SPECIAL POSITION OPERATORS
*
        IF(NSPE.NE.0) THEN
          WRITE(NWRITE,1095)(((INTEGER(SP(I,J,K),0),J=1,3),TSP(I,K),
     *          I=1,3), K , MUSP(K),K=1,NSPE)
 1095     FORMAT((1X,'SPOS',5X,2(3I2,F10.7,4X),3I2,F10.7,I3,I2))
          ENDIF
        WRITE(NWRITE,1050)
        ENDIF
*
* WRITE 'END' CARD
*
      WRITE(NWRITE,2000)
 2000 FORMAT('END',//)
      RETURN
      END 
c
c
c
c
      SUBROUTINE INCREA
C
C   ****SUBROUTINE****     INCREA
C                           ------
C WRITTEN 86/06/05 BY NINA HEINIG
C
C WRITES OUT SYMMETRY TABLE
C
      IMPLICIT CHARACTER*8 (G)
C
      COMMON /SYM/ NSYM,MSYM,S(3,3,48),T(3,48),NLT,TU(3,4),GSYM
      COMMON /FILES/ NIN,NOUT,NDEBUG,NREAD,NWRITE,NSPGP
      COMMON /GSYM/ GLAT,GCENT
C
      IF(NDEBUG .NE. 0) WRITE(NDEBUG, 2000)NSYM,NLT,TU
 2000 FORMAT(' *DEBUG**INCREA*2000 NSYM,NLT,TU:',2I4,4(/,1X,3F5.0)) 
      DO 10 K=1,NSYM
        WRITE(NWRITE,1000)((INTEGER(S(I,J,K),0),J=1,3),T(I,K),I=1,3),K
 1000   FORMAT(1X,'SYOP',5X,2(3I2,F10.7,4X),3I2,F10.7,I3)
   10   CONTINUE
      RETURN
      END
c
c
c
c
      subroutine conxyz(sp,tp,hbuf)
c
c     written by D. Altermatt
c     McMaster University, Nov. 1984
c
c     calls integer
c
c     for a given operator (sp/tp) this subroutine constructs the
c     general representation with x,y,z as in the International 
c     Tables for Crystallograpy and as required by DLS
c
      implicit character*80 (h)
      character*8 g1(12),g2(12),g3(12),g4(12)
c
      dimension sp(3,3),tp(3)
      data g1/'1/8','1/4','3/8','1/2','5/8','3/4','7/8','1',4*'   '/
      data g2/'1/12','1/6','1/4','1/3','5/12','1/2','7/12','2/3',
     1        '3/4','5/6','11/12','1'/
      data g3/'9/8','5/4','11/8','3/2','13/8','7/4','15/8','2',4*' '/
      data g4/'13/12','7/6','5/4','4/3','17/12','3/2','19/12','5/3',
     1        '7/4','11/6','23/12','2'/
c
      hline=' '
      do 10 i=1,3
        j=(i-1)*17+1
        if((abs(sp(i,1))+abs(sp(i,2))+abs(sp(i,3))+abs(tp(i))).lt.
     1    0.01) then
c   --    i-th coordinate 0
          hline(j:j)='0'
          if(i.ne.3) hline(j+1:j+1)=','
          goto 10
          endif
        if(integer(sp(i,1),0).eq.0) goto 20
c   --    X-component
          hline(j+1:j+1)='X'
          if(sp(i,1).gt.1.01) hline(j:j)='2'
          if(sp(i,1).lt.-0.1) hline(j:j)='-'
20      if(integer(sp(i,2),0).eq.0) goto 30
c   --    Y-component
          hline(j+4:j+4)='Y'
          if(sp(i,2).gt.1.01) hline(j+3:j+3)='2'
          if(sp(i,2).lt.-0.1) hline(j+2:j+2)='-'
          if(hline(j+1:j+2).eq.'X ') hline(j+2:j+2)='+'
30      if(integer(sp(i,3),0).eq.0) goto 40
c   --    Z-component
          hline(j+7:j+7)='Z'
          if(sp(i,3).gt.1.01) hline(j+6:j+6)='2'
          if(sp(i,3).lt.-0.1) hline(j+5:j+5)='-'
          if(hline(j+5:j+5).eq.' '.and.hline(j:j+4).ne.'     ')
     1      hline(j+5:j+5)='+'
40      if(integer(tp(i),3).eq.0) goto 50
c   --    translation
          if(hline(j:j+7).ne.'       ') hline(j+8:j+8)='+'
          if(tp(i).lt.-0.1) hline(j+8:j+8)='-'
          do 45 k=1,12
            xk=k
            k1=k*125
            k2=integer((xk*83.33),0)
            if(iabs(integer(tp(i),3)).eq.k1) 
     1        hline(j+9:j+14)=g1(k)(1:6)
            if(iabs(integer(tp(i),3)).eq.(k1+1000))
     1        hline(j+9:j+14)=g3(k)(1:6)
            if(iabs(integer(tp(i),3)).eq.k2) 
     1        hline(j+9:j+14)=g2(k)(1:6)
            if(iabs(integer(tp(i),3)).eq.(k2+1000))
     1        hline(j+9:j+14)=g4(k)(1:6)
45          continue
50        if(i.ne.3) hline(j+15:j+15)=','
10      continue
c
c     remove blanks, store in hbuf
c
      hbuf=' '
      j=1
      do 100 i=1,80
        if(hline(i:i).eq.' ') goto 100
        hbuf(j:j)=hline(i:i)
        j=j+1
100     continue
      return
      end
C
C
C
      SUBROUTINE MMINM(C,A,B,N)
C     -------------------------
C
C   WRITTEN 84.1.25 BY I.D.BROWN
C
C   CALCULATES C(N,N)=A(N,N)-B(N,N)
C
      DIMENSION A(N,N),B(N,N),C(N,N)
      DO 50 I=1,N
        DO 50 J=1,N 
   50     C(I,J)=A(I,J)-B(I,J)
      RETURN
      END 
*
*
*
*
      SUBROUTINE MXM (C,A,B,N)
*---- CALCULATES THE MATRIX PRODUCT FOR N*N MATRICES AS
*     A * B = C
*
*     A, B   INPUT MATRICES
*     C      OUTPUT MATRIX
*     N      ORDER OF MATRICES
*
*     WRITTEN BY D. ALTERMATT 
*---- MCMASTER UNIVERSITY, DAWN OF ORWELL'S YEAR
      DIMENSION A(N,N),B(N,N),C(N,N)
      DO 10 I=1,N
      DO 20 K=1,N
      C(I,K)=0
      DO 30 J=1,N
      C(I,K)=C(I,K)+A(I,J)*B(J,K)
   30 CONTINUE
   20 CONTINUE
   10 CONTINUE
      RETURN
      END 
C
C
C
C
      SUBROUTINE MXV(VPROD,SM,V)
C
C  *****SUBROUTINE*****     MXV
C                           ---
C
C   WRITTEN 83.9.7 BY I D BROWN
C
C   CALCULATES VPROD(3) = SM(3,3)*V(3)
C
      DIMENSION VPROD(3),SM(3,3),V(3)
      DO 50 I=1,3
         VPROD(I) = 0
         DO 50 J=1,3
   50       VPROD(I) = VPROD(I) + SM(I,J)*V(J)
      RETURN
      END 
*
*
*
*
*
      FUNCTION ICOMMA (A,B,N,ID)
*---- COMPARES TWO N*N MATRICES, GIVES ICOMMAT=0 IF THEY
*     ARE EQUIVALENT WITHIN ID DECIMAL PLACES, ELSE ICOMMAT=1
*
*     CALLS INTEGER
*
*     A, B  MATRICES TO CHECK 
*     N     ORDER OF MATRICES 
*     ID    NUMBER OF DECIMAL PLACES TO CHECK
*
*     WRITTEN BY D. ALTERMATT 
*     MODIFIED 84.5.16 BY I D BROWN
*---- MCMASTER UNIVERSITY, DAWN OF ORWELL'S YEAR
      DIMENSION A(N,N),B(N,N) 
      ICOMMA=0
      DO 10 I=1,N
      DO 20 J=1,N
        IF(INTEGER(A(I,J),ID).NE.INTEGER(B(I,J),ID)) THEN
          ICOMMA=1 
          RETURN
          ENDIF
   20 CONTINUE
   10 CONTINUE
   99 RETURN
      END 
C
C
C
C
      FUNCTION ICOMVE (X,Y,N,ID)
C     --------------------------
C
C   WRITTEN 84.2.7 BY I D BROWN
C
C   RETURNS THE VALUE OF 0 IF X(N) AND Y(N) ARE IDENTICAL TO ID DECIMAL
C   PLACES
C
      DIMENSION X(N),Y(N)
      ICOMVE=0
      DO 100 I=1,N
  100   IF(INTEGER(X(I),ID).NE.INTEGER(Y(I),ID)) ICOMVE=1
      RETURN
      END 
C
C
C
C
      FUNCTION INTEGER(X,N)
C
C   *****FUNCTION*****     INTEGER
C                          -------
C
C   WRITTEN 81.1.7 BY I D BROWN
C
C   RETURNS THE VALUE OF X*10**N AS THE NEAREST INTEGER
C
      Y = 1
      IF(N.LT.1) GO TO 20
      DO 10 I=1,N
   10    Y = Y*10.
   20 IF(X.LT.0.) GO TO 30
      INTEGER = X*Y + 0.5
      RETURN
   30 INTEGER = X*Y - 0.5
      RETURN
C**  THIS PROGRAM VALID ON FTN4 AND FTN5 **
      END 
C
C
C
*
      function second (itype)
*
*---  system dependent runtime call function, VAX-version
*
*     itype = 1      initialize system timer
*           = 2      print elapsed CPU time to SYS$OUTPUT
*           = 3      return elapsed CPU time
*
*     written by Daniel Altermatt
*---  McMaster University, September 1984
*
      external LIB$INIT_TIMER
      external LIB$STAT_TIMER
      external LIB$SHOW_TIMER
      character*80 htime
*
      second = 0.0
      if (itype.ne.1) goto 10
      call LIB$INIT_TIMER(it)
      return
*
10    if (itype.ne.2) goto 20
      icode = 2
      call LIB$SHOW_TIMER(it,icode,,htime)
      return
*
20    icode = 2
      call LIB$STAT_TIMER(icode,itime,it)
      second = itime / 100.
      return
      end
C
C
C
C
*     function second (itype)
C
C---  system dependent runtime call function, CYBER-version
C
C     itype = 1      initialize system timer
C           = 2      print elapsed CPU time to SYS$OUTPUT
C           = 3      return elapsed CPU time
C
C     written by Daniel Altermatt
C---  McMaster University, June 1985
C
*     EXTERNAL SECOND()
C
*     IF (ITYPE.LE.1) THEN
*        SECOND = 0.0
*     ELSE IF (ITYPE.EQ.2) THEN
*        WRITE(OUTPUT,1000) SECOND()
*1000    FORMAT(10X,'ELAPSED CPU-TIME SINCE START: ',F10.2,' SECONDS')
*     ELSE
*        SECOND = SECOND()
*        ENDIF
*     RETURN
*     END
C
C
C
C
      SUBROUTINE VMINV(Z,X,Y)
C
C   *****SUBROUTINE***** VMINV
C                        -----
C
C   WRITTEN 83.9.23 BY I D BROWN
C
C   CLACULATES Z(3) = X(3)-Y(3)
C
      DIMENSION X(3),Y(3),Z(3)
      DO 10 I=1,3
   10   Z(I) = X(I) - Y(I)
      RETURN
      END 
C
C
C
C
      SUBROUTINE VPLUSV(VPROD,V1,V2)
C
C   ******SUBROUTINE******    VPLUSV
C                             ------
C
C   WRITTEN 89.9.7 BY I D BROWN
C
C   CALCULATES VPROD(3) = V1(3) + V2(3) 
C
      DIMENSION VPROD(3),V1(3),V2(3)
      DO 50 I=1,3
   50    VPROD(I) = V1(I) + V2(I)
      RETURN
      END 


