     		PROGRAM TABLES
C
C . . THIS PROGRAM, IN CONJUNCTION WITH  THE CMD-PDS-GRAMPS
C . . ENVIRONMENT IN ABBOTT LABS,  PROVIDES THE CAPABILITY
C . . OF DISPLAYING AND/OR INTERACTING WITH THE SPACE GROUP SYMMETRY
C . . OF A PARTICULAR SPACE GROUP. IT CAN BE USED TO STUDY CRYSTAL PACKING,
C . . CRYSTAL CONTACTS,  AND AS A GUIDE IN UNDERSTANDING CRYSTALLOGRAPHIC
C . . SYMMETRY TO THE UNINITIATED, AS WELL AS THE PROFESIONAL CRYSTALLOGRAPHER.
C . .
C . . IT IS THE PURPOSE OF THIS PROGRAM TO COMPLEMENT PICTORIALLY  THE
C . .   'INTERNATIONAL TABLES FOR CRYSTALLOGRAPHY VOL. I'
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
C . .                                                                   .
C . .  . . . . . . . . . THIS IS VERSION  2.0  . . . . . . . . . . . . ..
C . .                     ABBOTT LABORATORIES 
C . .                     ABBOTT PARK, ILL., 60064                      . 
C . .                     C. ABAD-ZAPATERO
C . .                     PHONE (312)-937-0294
C . .  .  .  .  .  . . . . . . . . . . . . . . . . . . . . . . . . . . ..
C . . PROGRAM TABLES  VERSION 2.0  MARCH.11.1987.                       .
C . . VERSION THAT WILL BE RELEASED TOGETHER WITH THE J. APPLY. CRYSTALL.
C . . PAPER BY C. ABAD-ZAPATERO AND T. J. O'DONNELL..
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C     Modified 3/11/87  by TJ O'Donnell to produce GRAMPS output files,
C     making CMD/PDS (which only Abbott has) unnecessary.
C     GRAMPS output files are produced for the unit cell, symmetry operators
C     and the packing for the corresponding space group.
C . .
C . .
C . . VERSION TESTED UP TO HEXAGONAL SPACE GROUPS WITH SYMMETRY        .
C . . ELEMENTS AT 1/12,2/12  ETCC. '12' IN DENOMINATOR SHOULD BE
C . . REPLACED BY '9'.                                 .
C . .
C . . DIFFERENT SYMMETRY ELEMENTS (2-FOLDS,31-SCREWS,ETCC. ARE DEFINED  .
C . . BY USING DIFFERENT PSEUDO-ATOM TYPES IN THE CELL, AND OPS FILES    .
C . . ONE CAN USE COLOR BY TYPE AND HAVE THE DIFFERENT SYMMETRY OPERATORS
C . . WITH DIFFERENT COLORS:
C . .
C . . HEXAGONAL (SCREW) AXES 6,61,62,63 ETC.                  : BLUE(N) 
C . . TRIGONAL (SCREW) AXES 3,31,32 ETC.                      : YELLOW (S) 
C . . TETRAGONAL (SCREW) AXES 4,41,42 ETCC                    : BLUE (N)
C . . TWO FOLD (SCREW) AXES 2,21                              : RED  (O)
C . .
C . . SCREW AXES IN GENERAL WILL APPEARED PARTITIONED IF THE DATA
C . . FILES CONTAINING THE ORIGIN AND END OF THE VECTOR ARE SPECIFIED
C . . SO THAT ONLY ONE OF THE POINTS HAS THE CORRESPONDING N-FOLD LABEL. 
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C . .
C PROGRAM TO GENERATE A FULL CRYSTAL CELL USING THE CRYSTALLOGRAPHIC
C SYMMETRY OPERATORS OPERATING ON A MOTIF FILE
C COMING FROM INTERACT OR FROM GRAMPS 
C OR CREATED BY HAND FOR THE PURPOSE OF CHECKING THE PACKING
C AND/OR THE SYMMETRY OF THE SPACE GROUP
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C ..........
C . . . . .INPUT FILES:
C . .
C . LUN=1  SPACE GROUPS ARE INPUT AS STRINGS (MAXIMUM 36 CHARACTERS)
C . . . . "A LA INTERNATIONAL TABLES " : X,Y,Z;-X,Y,-Z ETCC. TO
C . . . .  A TOTAL OF NEQ STRINGS, ONE PER LINE. THE FILE ALSO
C . . . .  CONTAINS VECTORS DEFINING THE POSITION OF THE DIFFERENT
C . . . .  SYMMETRY OPERATORS IN THREE DIMENSIONS. THESE VECTORS ARE
C . . . .  SPECIFIED IN A FORMAT COMPLETELY ANALOGOUS TO THE ONE
C . . . .  USED IN THE INTERNATIONAL TABLES
C . . . .
C . LUN=2 IS MOTIF FILE IN ABBOTT'S ***.MSF FORMAT 
C . . . .  LABEL,X,Y,Z,W,ATOMTYPE
C . . . .  CONTAINING THE MOTIF. THIS MOTIF COULD BE A MOLECULE
C . . . .  OR A WORDLY OBJECT, AS LONG AS IT HAS THE .MSF FORMAT
C . . . .  (SEE ROUTINES HEADMSF (READS HEADER MSF FORMAT) AND
C . . . .   OUTCELL, OUTOPS
C . . . . 
C . .
C . . LUN=3 MOTIF.FILE IN BROOKHAVEN DATA BANK FORMAT. IN GENERAL
C . . . . . ONLY CA*S COORDINATES WILL BE EXTRACTED FOR THIS
C . . . . . TYPE OF FILE SO THAT THE NUMBER OF ATOMS IN THE
C . . . . . SCREEN IS MANAGEABLE.  
C . . .
C . . LUN=4 MOTIF FILE IN CAMBRIDGE DATA BANK FORMAT.
C . . . . . .
C . . . OUTPUT FILES
C . .
C . . LUN=7 IS THE (SCRATCH) OUTPUT FILE IN CAMBRIDGE DATA BANK FORMAT
C . . RECORD 1:  TITLE (80A)
C . . RECORD 2:  CELL CONSTANTS A,B,C,ALFA,BETA,GAMMA
C . . RECORD 3:  SYMMETRY OPERATORS (THERE ARE NEQ (NUMBER OF EQUIV. POS))
C . . RECORD 4,ETCC... COORDINATES OF ATOMS IN THE A.U. IN FRACTIONAL       
C . . COORDINATES IN A LATTICE CONSISTENT WITH THE  
C . . THE ONE TO BE USED IN THE FILE IN GRID FILE (LUN=8).
C . .
C . . LUN=8  UNIT CELL PARALLELEPIPED WITH PSEUDO-ATOMS AS VERTICES
C . .        IN .PIX FORMAT READY TO BE DISPLAYED BY GRAMPS
C . . LUN=9  VECTORS DEFINING LOCATION OF SYMMETRY OPERATORS
C . .        IN .PIX FORMAT READY TO BE DISPLAYED BY GRAMPS
C . . LUN=10 CRYSTAL FILE CONTAINING THE MOTIF IN LUN=1 AND THE
C . .        CORRESPONDING ONES GENERATED BY THE CRYSTAL SYMMETRY.
C . .        THIS FILE NEEDS TO BE CONVERTED TO .PIX FORMAT TO DISPLAY IT.
C . . 
C . . . . . . . . .
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      PARAMETER MAXATOM=2000
      PARAMETER MAXSYM =32
      PARAMETER MAXOPS =200
C . .
C . .
C . .
      LOGICAL ALLBUTH,ONLYCA,OKAY
      CHARACTER*64 INFIL1,INFIL2,OUTFIL1,OUTFIL2,OUTFIL3,OUTFIL4
      CHARACTER*80 LIN1,LIN2,LATLAB
      INTEGER*2   IBONDS
      CHARACTER*4 BONDS
      CHARACTER*6 IATN
      CHARACTER*8 DATE,TIME
      CHARACTER*9 PROG
      CHARACTER*9 VERS
      CHARACTER*2 CA,SYMEL,BLANK2
      CHARACTER*1 IATO,H,N,C,O,S,LT,BLANK1,ANSWER,IATID
      CHARACTER*40 SPG,FULSYM 
      DATA DATE,TIME  /'11-03-87','14:30:00'/
      DATA PROG,VERS  /' TABLES  ','VERS 2.0 '/
      DATA H,N,C,O,S,BLANK1 /'H','N','C','O','S',' '/
      DATA BLANK2 /'  '/
C . . DIMENSION VECTORS
      DIMENSION TRS(3),VI(3),VO(3),RF(3),RO(3),XOR(3,4)
C . . DIMENSION OPERATORS
C . . CEL1(6) REAL CELL CONSTANTS A,B,C,ALFA,BETA,GAMMA (COSINES)
      DIMENSION CEL1(6),CEL2(6),SYM(3,4,MAXSYM),SYOP(12,MAXSYM)
      DIMENSION TM(3,3),TMI(3,3),UNI(3,3),R(3,4),WORK(3,4)
C . . V000 CONTAINS THE ORTHOGONAL ANGSTROMS COORDINATES OF THE VERTICES
C . . OF THE UNIT CELL (IN COMMON BLOCK).
C . . G000 CONTAINS THE FRACTIONAL COORDINATES FOR
C . . THE VERTICES          
C . . BONDS CONTAINS THE INFORMATION TO JOIN THE VERTICES OF THE UNIT CEL
C . . (IN COMMON BLOCK).
      DIMENSION G000(3,8)
C . .
C . . SET COMMON BLOCK FOR THE CELL FILE:
C . . 
      COMMON /CELL/V000(3,8),BONDS(2,12),IBONDS(2,12),C,WG1,BONTY1     
C . .
C . . SET COMMON BLOCK FOR THE SYMOPS FILE : MAXOPS = 200
      COMMON /SYMOPS/SVEC(3,MAXOPS),IATID(MAXOPS),WG2,BONTY2
C . .
      LOGICAL*1 SYMOP(36)
      DATA OPT1/'TITL'/
      DATA OPT2/'CELL'/
      DATA OPT3/'SYMM'/
C . . BOND TYPES WILL BE SET UP AS ONE (1.0; i.e. single bonds)
C . . FOR THE CELL, AND SYMMETRY OPERATORS..
      DATA BONTY1,BONTY2/1.0,1.0/ 
C . . VERTICES OF UNIT CELL IN FRACTIONAL COORDINATES
      DATA G000(1,1)/0./,G000(2,1)/0./,G000(3,1)/0./
      DATA G000(1,2)/1./,G000(2,2)/0./,G000(3,2)/0./
      DATA G000(1,3)/0./,G000(2,3)/1./,G000(3,3)/0./
      DATA G000(1,4)/1./,G000(2,4)/1./,G000(3,4)/0./
      DATA G000(1,5)/0./,G000(2,5)/0./,G000(3,5)/1./
      DATA G000(1,6)/1./,G000(2,6)/0./,G000(3,6)/1./
      DATA G000(1,7)/0./,G000(2,7)/1./,G000(3,7)/1./
      DATA G000(1,8)/1./,G000(2,8)/1./,G000(3,8)/1./
C . . BONDS DESCRIPTION FOR UNIT CELL
      DATA BONDS(1,1)/'C100'/,BONDS(2,1)/'C101'/,
     .     BONDS(1,2)/'C100'/,BONDS(2,2)/'C102'/,
     .     BONDS(1,3)/'C101'/,BONDS(2,3)/'C103'/,
     .     BONDS(1,4)/'C102'/,BONDS(2,4)/'C103'/,
     .     BONDS(1,5)/'C104'/,BONDS(2,5)/'C105'/,
     .     BONDS(1,6)/'C104'/,BONDS(2,6)/'C106'/,
     .     BONDS(1,7)/'C105'/,BONDS(2,7)/'C107'/,
     .     BONDS(1,8)/'C106'/,BONDS(2,8)/'C107'/,
     .     BONDS(1,9)/'C100'/,BONDS(2,9)/'C104'/,
     .     BONDS(1,10)/'C101'/,BONDS(2,10)/'C105'/,
     .     BONDS(1,11)/'C102'/,BONDS(2,11)/'C106'/,
     .     BONDS(1,12)/'C103'/,BONDS(2,12)/'C107'/
C . .  BONDS DESCRIPTION FOR PIX FILE FOR UNIT CELL
      DATA IBONDS(1,1)/0/,IBONDS(2,1)/1/,
     .     IBONDS(1,2)/0/,IBONDS(2,2)/2/,
     .     IBONDS(1,3)/1/,IBONDS(2,3)/3/,
     .     IBONDS(1,4)/2/,IBONDS(2,4)/3/,
     .     IBONDS(1,5)/4/,IBONDS(2,5)/5/,
     .     IBONDS(1,6)/4/,IBONDS(2,6)/6/,
     .     IBONDS(1,7)/5/,IBONDS(2,7)/7/,
     .     IBONDS(1,8)/6/,IBONDS(2,8)/7/,
     .     IBONDS(1,9)/0/,IBONDS(2,9)/4/,
     .     IBONDS(1,10)/1/,IBONDS(2,10)/5/,
     .     IBONDS(1,11)/2/,IBONDS(2,11)/6/,
     .     IBONDS(1,12)/3/,IBONDS(2,12)/7/
CC
C . .
C . .
C . . T-MATRIX TRANSFORMS FROM FRACTIONAL OBLIQUE COORDINATES
C . . TO ANGSTGROMS IN ORTHOGONAL AXES
C . .
C . . MATRIX T- INVERS DOES THE OPPOSITE: ORTHOGONAL ANGSTROMS
C . . TO FRACTIONAL OBLIQUE
C . . A CHECK IS DONE FOR TM*TMI = UNITY MATRIX
C . .
C . . DEFINE I/O FILES
C . .
C . . INPUT FILE CONTAINING SPACE GROUP OPERATORS "A LA INTERNATIONAL
C . .  TABLES": X,Y,Z; -X,Y,-Z ETC. AS MANY AS NEQ STRINGS ONE/LINE
C . . TO GENERATE THE CORRESPONDING SYMMETRY
C . . FILE [ABAD.INTL.TABLES] SPG.DAT 
       INTL=1
C . .
C . . 
C . . INPUT FILE CONTAINING MOTIF IN .MSF FORMAT
       IMSF = 2
C . . INPUT FILE MOTIF IN BROOKHAVEN DATA BANK 
       IPDB = 3
C . . INPUT FILE MOTIF IN CAMBRIDGE DATA BANK (NOT WORKING YET. DEC.26.1986)
       ICAM = 4
C . .
C . .
C . . OUTPUT FILE CONTAINING MOTIF AND SYMMETRY OPERATORS  LUN=2 .CRY FORMAT
C . . IN FRACTIONAL COORDINATES TO BE USED BY CALL XTAL (OPTIONAL)
C . . TO GENERATE THE CRYSTAL PACKING (SEE BELOW).
        JCRY=7
C . . OUTPUT FILE CONTAINING GRID POINTS TO FORM UNIT CELL(S) .MSF FORMAT
        IGRD=8
C . . OUTPUT FILE CONTAINING VECTORS DEFINING THE SYMMETRY ELEMENTS IN THE
C . . UNIT CELL IN ORTHOGONAL COORDINATES .MSF FORMAT READY FOR INTERACT
        ISO=9
C . .OUTPUT FILE CONTAINING THE CREATED CRYSTAL PACKING IN CAMBRIDGE FORMAT
C . .AND READY FOR CONVERT ====> INTERACT ==== DISPLAY
        ICRY=10
C . .
C . . . SET NUMBER OF VERTICES AND NUMBER OF BONDS TO DEFINE UNIT CELL
        NVT = 8
        NBD = 12
C . . ZERO OUT SYMMETRY OPERATORS ARRAY
        DO I=1,3
           DO J=1,4
            DO K=1,MAXSYM
              SYM(I,J,K) =0.0
            END DO
           END DO
        END DO
C . . ZERO OUT ORIGINS ARRAY FOR DIFFERENT CELL TYPES 
         DO I=1,3
          DO J=1,4
           XOR(I,J) =0.0
          END DO
         END DO
C . . ZERO OUT SYMMETRY VECTORS ARRAYS: 
         DO I=1,3
          DO J=1,MAXOPS
           SVEC(I,J) =0.0
          END DO
         END DO
		
C. . JUST A FEW CONSTANTS 
        D2R = 3.14159265/180.
        R2D = 1./D2R
        WG1 = 0.1000
        WG2 = 0.1000

C . . . .
	TYPE 9999
 9999   FORMAT(5X,'--------------------------------------------'/,
     .  5X,'-------- INTERNATIONAL TABLES ---------------------'/,
     .  5X,'       STEREO GENERATION PROGRAM VERSION 2.0       '/,
     .  5X,'   MAR.11.1987. BY C. Abad-Zapatero ABBOTT LABS    '/,
     .  5X,'   MAY.12.1987. BY TJ O''Donnell     ABBOTT LABS    '/,
     .  5X,'    CORRESPONDING TO THE J. APPL. CRYST. PAPER     '/,
     .  5X,'---------------------------------------------------'/)   
C . . . SPECIFY INPUT AND OUTPUT FILES   
C . .
        TYPE 1
1       FORMAT(1X,'ENTER NAME OF INPUT FILE1: [SPG.DAT] ')
        ACCEPT 2,INFIL1
2       FORMAT(A)
        TYPE 4,INFIL1
4       FORMAT(' INPUT FILE 1: [SPG.DAT]',A/' ENTER NAME FOR ',
     . 'OUTPUT FILE 1: [SPG1.CRY] ')
        ACCEPT 2, OUTFIL1
        TYPE 5, OUTFIL1   
5       FORMAT(' OUTPUT FILE 1: SPG1.CRY ',A/' ENTER NAME FOR',
     . ' OUTPUT FILE 2 : [CELL.PIX] ')
        ACCEPT 2,OUTFIL2
        TYPE 6,OUTFIL2
6       FORMAT(' OUTPUT FILE 2: CELL.PIX ',A/' ENTER NAME FOR',
     . ' OUTPUT FILE 3 : [SYMOPS.PIX] ')
        ACCEPT 2,OUTFIL3
        TYPE 7,OUTFIL3
7       FORMAT(' OUTPUT FILE 3: [SYMOPS.PIX] ',A)
C        TYPE 8
C8       FORMAT(' LABEL OF THE FIRST RECORD OF CELL.MSF FILE? ')
C        ACCEPT 2, LATLAB
C . . . . OPEN FILES . . 
C . . INPUT
        OPEN(UNIT=INTL,NAME=INFIL1,STATUS= 'OLD',READONLY
     .,DEFAULTFILE='SPG.DAT')
C . . OUTPUT
        OPEN(UNIT=JCRY,NAME=OUTFIL1,STATUS='NEW',SHARED
     .,DEFAULTFILE='SPG1.CRY')
        OPEN(UNIT=IGRD,NAME=OUTFIL2,STATUS='NEW'
     .,DEFAULTFILE='CELL.PIX')
        OPEN(UNIT=ISO,NAME=OUTFIL3,STATUS='NEW'
     .,DEFAULTFILE='SYMOPS.PIX')
C . . . .        
C . . . . UNIT CELL DATA FOR ORTHOGONALIZATION AND/OR
C . . . . DEORTHOGONALIZATION
        WRITE(6,6000) 
 6000   FORMAT(5X, 'PLEASE TYPE UNIT CELL PARAMETERS  '/)
        READ(5,* )  CEL2
         WRITE(6,200) CEL2
 200     FORMAT(5X, 'UNIT CELL CONSTANTS FOR ORTHO/DEORTHO ',/,
     .5X, 'A , B, C, ALF, BET, GAM ', 6F10.5)
            DO K=1,6
               CEL1(K) = CEL2(K)
            END DO
	    CEL1(4) = COS(D2R*CEL1(4))
	    CEL1(5) = COS(D2R*CEL1(5))
	    CEL1(6) = COS(D2R*CEL1(6))
C . . COMPUTE VOLUME
            VOL = VOLUME(CEL1)
C. . .
C . . . 
C . . 
	    WRITE(6,10) VOL
 10     FORMAT(/,5X, 'VOLUME OF UNIT CELL ',/,10X,F12.2,1X,'A**3',/)
         WRITE(6,210) 
 210     FORMAT(5X, 'DE-ORTHOGONALIZATION DONE SO THAT 
     .  THE X- AXIS REMAINS UNCHANGED ',/,
     ./,5X, 'DE-ORTHO MATRIX ', /)
        CALL TMATRIX(CEL1,VOL,TM)
        DO 65 K=1,3
        WRITE(6,37) (TM(K,L),L=1,3)
  65    CONTINUE
  37    FORMAT(/,5X, 'T-MATRIX   ', 3F12.5)
        CALL INVERS (TM,TMI)
        DO 75 K=1,3
        WRITE(6,77) (TMI(K,L),L=1,3)
  75    CONTINUE
  77    FORMAT(/,5X, 'T-MATRIX INVERS', 3F12.5)
C. . CHECK OF INVERS TRANSFORMATION
C. . MATRIX PRODUCT OF TM*TMI
        CALL MATPRO(UNI,TM,TMI)
C . . . . . . . . .. . . .    
        DO 223 K=1,3
        WRITE(6,225) (UNI(K,L),L=1,3)
225     FORMAT(/,5X, 'TM * TMI = ', 3F8.4,/)
223     CONTINUE
C . . . .
        WRITE(6,1030)INFIL2 
 1030   FORMAT(5X,'READING SPACE GROUP DATA ',A/)
        READ(INTL,3000) LIN1
        WRITE(6,1037) LIN1
        READ(INTL,1040)NSPG,NEQ,SPG
        DO K=1,10
         IF(SPG(K:K).NE.BLANK1) GO TO 4880
        END DO
 4880   LT = SPG(K:K)
C . . SET ILT VARIABLE FOR FUTURE USE IN SETTING UP NOR AND XOR(3,NOR)
        IF (LT.EQ. 'P' ) ILT=1
        IF (LT.EQ. 'C' ) ILT=2
        IF (LT.EQ. 'F' ) ILT=3
        IF (LT.EQ. 'I' ) ILT=4
        IF (LT.EQ. 'A' ) ILT=5
C . .
C . .
        GO TO (4300,4310,4320,4330,4340), ILT
C . .   PRIMITIVE LATTICE
 4300   NOR = 1
C . .   XOR ALREADY ZEROED OUT
        GO TO  4350
C . .   C-CENTERED
 4310   NOR = 2
        XOR(1,2 )= 0.5
        XOR(2,2 )= 0.5
        XOR(3,2 )= 0.0
        GO TO 4350
C . .   F-CENTERED
 4320   NOR=4
        XOR(1,2)=0.
        XOR(2,2)=0.5
        XOR(3,2)=0.5
        XOR(1,3)=0.5
        XOR(2,3)=0.0
        XOR(3,3)=0.5
        XOR(1,4)=0.5
        XOR(2,4)=0.5
        XOR(3,4)=0.
        GO TO 4350
C . .   I-CENTERED
 4330   NOR=2
        XOR(1,2) =0.5
        XOR(2,2) =0.5
        XOR(3,2) =0.5
        GO TO 4350
C . .   A-CENTERED LATTICE
 4340   NOR=2
        XOR(1,2) =0.
        XOR(2,2) = 0.5
        XOR(3,2) = 0.5
        GO TO 4350
C . . END OF SETTING UP NUMBER OF ORIGINS AND ORIGIN ARRAYS
 4350   CONTINUE      
        IF (NEQ.LE.0 ) GO TO 5000
        WRITE(6,1050) SPG,NEQ,LT,ILT
C . .
C . . TOTAL NUMBER OF EQUIVALENT POSITIONS NEQT = NEQ * NOR
        NEQT = NEQ*NOR
        IF( NEQ.GT.0.AND.NEQT.LT.MAXSYM) GO TO 5000
C . .
        WRITE(6,1060)
        STOP 'NUMBER OF EQUIVALENT POSITIONS GT MAXSYM'
 5000   CONTINUE
        IF (NEQ.GT.0) GO TO 5005
C . .
        STOP 'BEWARE: NEQ .LE. 0 IN YOUR SPG.DAT FILE'
C. . .
C. . . START TO BUILD THE SPACE SPG1.CRY FILE
C. . 
 5005   WRITE(JCRY,1055)OPT1, SPG
 1055   FORMAT(A,5X,10A4)
        WRITE(JCRY,1065)OPT2,CEL2
 1065   FORMAT(A,4X,6F8.3)
        IF (NEQ.EQ.0 ) GO TO 5010    
        WRITE(6,1035) 
 1035   FORMAT(5X,'READING SYMMETRY EQUIVALENT POSITIONS:
     .   X,Y,Z ETC..',/) 
          DO IS = 1, NEQ
           READ(INTL,1070) SYMOP
           J = MATSYM(SYM(1,1,IS),SYMOP,ICOL)
           IF (J.LE.0) GO TO 2150
           WRITE(6,1071) J,I,SYMOP,ICOL
           STOP 'PROBLEMS WITH EQUIVALENT POSITIONS'
 2150      WRITE(6,1072) IS,SYMOP,((SYM(J,K,IS),K=1,4),J=1,3)
          END DO
 1037     FORMAT(5X,'HEADER RECORD IN [SPG.DAT]: ',A)
          READ(INTL,1040) NSYP,NSYV,FULSYM
          DO IS = 1, NSYP
C . . . READ FROM INTERNATIONAL TABLES DATA LOCATION OF SYMMETRY ELEMENT
C . . .  SYMOP(36) AND TYPE SYMEL(2) DEFAULT FOR SYMEL IS BLANK 
          READ(INTL,1080) SYMOP,SYMEL
C . . . DEFAULT SYMMETRY ELEMENT IS BLANK THEN PSEUDO ATOMS FOR
C . . . COLOR CODE ARE CARBONS OR COLOR GREEN WHEN COLORED BY TYPE
            IF (SYMEL.EQ. BLANK2 ) IATID(IS) = C       !GREEN COLOR DEFAULT 
            IF (SYMEL.EQ. '2 '   ) IATID(IS) = O       !RED COLOR 2-FOLDS
            IF (SYMEL.EQ. '31'   ) IATID(IS) = S       !YELLOW COLOR 31-SCREWS
            IF (SYMEL.EQ. '3 '   ) IATID(IS) = S       !YELLOW COLOR 3-FOLDS
            IF (SYMEL.EQ. '32'   ) IATID(IS) = S       !YELLOW FOR 32-ALSO
            IF (SYMEL.EQ. '21'   ) IATID(IS) = O       !RED COLOR 21-SCREWS
            IF (SYMEL.EQ. '41'   ) IATID(IS) = N       !BLUE COLOR 4-FOLDS
            IF (SYMEL.EQ. '4 '   ) IATID(IS) = N       !41,42, OR 43 SCREWS
            IF (SYMEL.EQ. '42'   ) IATID(IS) = N
            IF (SYMEL.EQ. '43'   ) IATID(IS) = N      
            IF (SYMEL.EQ. '6 '   ) IATID(IS) = N       !BLUE COLOR 6-FOLDS 
            IF ((SYMEL.EQ. '61') .OR. (SYMEL.EQ.'63') 
     .  .OR. (SYMEL.EQ. '65'))     IATID(IS) = N 
            IF ((SYMEL.EQ. '62') .OR. (SYMEL.EQ.'64'))
     .                IATID(IS) = N 
          J = VECSYM(VI,SYMOP,ICOL)
          IF (J.LE.0) GO TO 2155
          WRITE(6,1171) J,I,SYMOP,ICOL
          STOP 'PROBLEMS WITH VECTORS SYMMETRY OPERATORS'
 2155     WRITE(6,1172) IS,SYMOP,SYMEL,(VI(K),K=1,3)
C . . . TRANSFORM FRACTIONAL COORDINATES OF VECTOR 1,0,0.333
C . . . TO ORTHOGONAL COORDINATES  FOR SYMOP.PIX FILE
        
          CALL MATVEC2(VO,TM,VI)
            DO LL=1,3
              SVEC(LL,IS) = VO(LL)
            END DO   
          END DO
          
 102      FORMAT(A,4X,3(3F4.0,1X,F7.5,1X))
 5010     CONTINUE
C . . SPACE GROUP SYMMETRY SPECIFIED BY S.G. NUMBER AND/OR SYMBOL
          WRITE(6,5020) NSPG,SPG,NEQ
 5020     FORMAT (5X, 'INTERNATIONAL TABLES S.G. NUMBER ', I5,/
     .           ,5X, 'INTERNATIONAL TABLES S.G. SYMBOL ', 10A4,/
     .           ,5X, 'NUMBER OF EQUIVALENT POSITIONS   ',  I5,/)  
C
C . . EXPAND THE NUMBER OF SYMMETRY OPERATORS TO INCLUDE THE 
C . . THE DIFFERENT ORIGINS FOR A-,C-,F- AND I- CENTERING (SEE ABOVE)
C . . FOR SETTING NOR (NUMBER OF ORIGINS) AND XOR(3,NOR) (ORIGIN COORDINATES)
C . .  
           DO  IS=1,NEQ
            DO IOR =1,NOR 
              DO K=1,3
                 DO L=1,3
                  SYM(K,L,IS+NEQ) = SYM (K,L,IS)
                 END DO 
                  SYM(K,4,IS+NEQ) = SYM (K,4,IS) + XOR (K,IOR)
              END DO
            END DO
           END DO
C
           DO IS=1,NEQT
           IP = 0
            DO K=1,3
               DO L=1,4
               IP = IP + 1
               SYOP(IP,IS) = SYM(K,L,IS)
              END DO
            END DO
           END DO
C . . WRITE OUT SYMMETRY OPERATORS ON SCRACTH FILE 
          DO L=1,NEQT
          WRITE(JCRY,1075) OPT3,( SYOP(K,L),K=1,12)
 1075     FORMAT(A,4X,3(3F4.0,1X,F7.5,1X))
          END DO
     
C
C . .
3000     FORMAT(A)
1040     FORMAT (2I5,A)
1050     FORMAT (5X,'SPACE GROUP ',A/5X,'EQUIVALENT POSITIONS',I5
     .  ,1X,'LATTICE TYPE ',1X,A,1X,' ILT ',I5)
1060     FORMAT (5X,'INVALID NUMBER OF POSITIONS')
1070     FORMAT (36A1)
1080     FORMAT (36A1,2X,A)
1071     FORMAT ('0*****MATSYM ERROR #'I3,' IN SYMMETRY OPERATOR'I2,
     .        '-- '36A1/' AT CHARACTER POSITION 'I2)
1072     FORMAT (//' SYMMETRY OPERATOR #'I3,5X,36A1,/
     .          3(5X,3F5.1,5X,F10.5/) )
1171     FORMAT ('0*****VECSYM ERROR #'I3,' IN SYMMETRY OPERATOR'I2,
     .        '-- '36A1/' AT CHARACTER POSITION 'I2)
1172     FORMAT (//' SYMMETRY VECTOR #'I3,5X,36A1,A/
     .          3(5X,3F10.5/) )
C ....
C . . . GENERATE AND STORE THE ORTHOGONAL COORDINATES OF THE 8-VERTEX
C . . . OF THE UNIT CELL BY MULTIPLYING  TM(3,3)*G000(3) 
C . . .
        DO K=1,NVT
           DO L=1,3
             VI(L) = G000(L,K)
           END DO
                 CALL MATVEC2(VO,TM,VI)
           DO M=1,3
             V000(M,K) = VO(M)
           END DO
C. . DEBUG
C       WRITE(6,7995) VI,TM,VO
C7995   FORMAT(1X,3(F6.2,1X),9(F8.2,1X),3(F6.2,1X))
C . . END DEBUG
         END DO
C . .
C . . WRITE CELL FILE 
C . .
C . . IGRD OUTPUT FILE NVT= NUMBER OF VERTICES 8
C . .                  NBD= NUMBER OF BONDS 12  FOR PARALLELEPIPED
C      CALL OUTCELL
C     . (IGRD,LATLAB,LIN2,DATE,TIME,PROG,VERS,IABOT,NVT,NBD)
      CALL OUTCELLPIX
     . (IGRD,DATE,TIME,PROG,VERS,IABOT,NVT,NBD)
C . .
C . . WRITE SYMOPS FILE 
C . .
C . . IGRD OUTPUT FILE NVT= NUMBER OF VERTICES 8
C . .                  NBD= NUMBER OF BONDS 12  FOR PARALLELEPIPED
C      CALL OUTOPS
C     . (ISO,LATLAB,LIN2,DATE,TIME,PROG,VERS,IABOT,NSYP,NSYV)
      CALL OUTOPSPIX
     . (ISO,DATE,TIME,PROG,VERS,IABOT,NSYP,NSYV)
C . .
C . . . .  
        CLOSE (UNIT=IGRD)
        CLOSE (UNIT=ISO)
C . . . CELL AND SYMOPS FILES COMPLETED.
C . . . 
C . . . DO YOU WANT TO GENERATE THE CRYSTAL SYMMETRY AND PACKING
C . . .
        TYPE 7000
 7000   FORMAT(/,5X,'UNIT CELL AND SYMMETRY OPERATORS DONE',/
     .  5X, 'DO YOU WANT TO GENERATE THE CRYSTAL? [Y/N]',/)
C . .
        READ(5,3000) ANSWER
         IF (ANSWER.EQ. 'Y'.OR. ANSWER.EQ. 'y' ) THEN

C . .  
        TYPE 2880
 2880   FORMAT(5X, 'SPECIFY THE FORMAT OF YOUR MOTIF ',/
     .,5X,' MSF TYPE FILE  MOTIF: 1',/
     .,5X,' BROOKHAVEN DATA BANK MOTIF: 2',/
     .,5X,' CAMBRIDGE DATA BANK MOTIF: 3 ',/
     .,5X,' PLEASE ENTER MOTIF TYPE===> ')
        READ(5,*) MOTIF
C . . BRANCH OFF ACCORDING TO MOTIF FILE TYPE:

       TYPE 3 
 3     FORMAT(5X, 'ENTER NAME OF YOUR MOTIF FILE: '/
     .,5X, ' MOTIF = 1   [MOTIF.MSF]  ',/
     .,5X, ' MOTIF = 2   [MOTIF.PDB]  ',/
     .,5X, ' MOTIF = 3   [MOTIF.CDB]  ')
       ACCEPT 2,INFIL2 
        TYPE 7010
 7010   FORMAT(5X ,'SPECIFY OUTPUT FILE NAME [SPG2.CRY]',/)
        READ(5,3000) OUTFIL4 
        WRITE(6,7015) 
 7015   FORMAT(/,5X, 'DO YOU WANT TO ROTATE AND/OR
     .  TRANSLATE YOUR MOTIF? [Y/N]',/)
        READ(5,3000) ANSWER
        IF (ANSWER.EQ.'Y'.OR.ANSWER.EQ.'y') THEN
          WRITE(6,6010)
 6010     FORMAT(5X,'  GIVE CAPA,PSI AND PHI ',/ 
     .  5X,'    PSI FROM Z-AXIS;PHI FROM:+X ',/)
	  READ(5,* ) CAPA,PSI,PHI
          WRITE(6,6020) 
 6020     FORMAT(5X,'TRANSLATION COMPONENTS? X,Y,Z (FRACT)',/)
        READ(5,*) TRS
C . .
C . . . .TRANSFORM FRACTIONAL COORDINATES OF TRANSLATION VECTOR
C . . . .INTO ORTHOGONAL COORDINATES
         ELSE
          CAPA = 0.0
          PSI = 0.
          PHI = 0.
          DO I0=1,3
           TRS(I0) = 0.0
          END DO
         END IF  
        
         CALL MATVEC2(VO,TM,TRS)
            DO LL=1,3
              TRS(LL) = VO(LL)
            END DO   
C .
C . . . .
        WRITE(6,250) CAPA,PSI,PHI,TRS
 250    FORMAT(5X, 'ROTATION ABOUT THE AXIS OF ', F6.2,/
     .,5X, 'AXIS DIRECTION DEFINED BY ', 'PSI = ',F10.2
     .,5X, 'PHI = ', F10.2, 2X, 'POLAR CO-ORDINATES ',/
     .,5X, 'TRANSLATION VECTOR X,Y,Z ANGSTROMS', 3F8.4,/ 
     .,5X, 'MATRIX OF THE TRANSFORMATION ',/)
C . . FORM THE ROTATION MATRIX FOR  CAPA DEGREES ABOUT AXIS
          CAP=D2R*CAPA
          PS =D2R*PSI
          PH =D2R*PHI
          CALL ROTMAT(CAP,PS,PH,R)
C . . .
          DO 290 K=1,3
          R(K,4) = TRS(K)
          WRITE(6,300) (R(K,L),L=1,4) 
 300      FORMAT(5X,   4F8.4,/)
 290      CONTINUE
       GO TO (7770,7775,7780), MOTIF

C . .   
C . .  READ HEADER RECORDS OF .MSF TYPE INPUT FILE BEFORE
C . .  READING THE X,Y,Z COORDINATES; OBTAIN
C . .  BASIC INFORMATION: NATOMS, NBONDS
C . .  AND OTHER: DATE,TIME,PROG,VERS,IABOT ETCC. 
C . .
C . . GET HEADER OF INPUT MSF TYPE FILE
C .
C . .
 7770    OPEN(UNIT=IMSF,NAME=INFIL2,STATUS='OLD',READONLY
     .,DEFAULTFILE='MOTIF.MSF')

       CALL HEADMSF(IMSF, LIN1,LIN2,DATE,TIME,PROG,VERS,IABOT,
     .             NUMATO,NUMBOND)
        WRITE(6,2995)
 2995   FORMAT(5X, ' ALL BUT HYDROGEN ATOMS OUTPUT ? Y/N  ',$)
        READ(5,3000) ANSWER
        IF (ANSWER.EQ. 'Y'.OR. ANSWER.EQ. 'y' ) THEN
          ALLBUTH = .TRUE.
        ELSE
          ALLBUTH = .FALSE.
          WRITE(6,2996)
 2996     FORMAT(5X, ' OXYGEN ATOMS WILL ALSO BE EXCLUDED  ')
        END IF
C . . NOW READ THE ATOMS AND WRITE THEM TO SCRATCH FILE JCRY
C . . IN FRACTIONAL COORDINATES FOR CRYSTAL GENERATION
C . .
        CALL ATOMSIN (ALLBUTH,IMSF,JCRY,NUMATO,R,TM,TMI,XOR,NOR)
C . .
        GO TO 7785
C . . AS BEFORE BUT FOR A PROTEIN DATA BANK MOTIF FILE
 7775   OPEN(UNIT=IPDB,NAME=INFIL2,STATUS='OLD',READONLY
     .  ,DEFAULTFILE='MOTIF.PDB')
        WRITE(6,8995)
 8995   FORMAT(5X, ' ONLY CA CARBONS? Y/N  ',$)
        READ(5,3000) ANSWER
        IF (ANSWER.EQ. 'Y'.OR. ANSWER.EQ. 'y' ) THEN
          ONLYCA = .TRUE.
        ELSE
          ONLYCA = .FALSE.
          WRITE(6,8996)
 8996     FORMAT(5X, 'PACKING USING ALL ATOMS: BEWARE!!')
        END IF
C . . READ PROTEIN DATA BANK FILE, EXAMINE HEADER RECORDS AND
C . . WRITE ATOMS INTO SCRATCH FILE JCRY FOR CRYSTAL GENERATION
C . . (ONLY CA CARBONS (DEFAULT) FOR PACKING STUDIES)
C . . IN FRACTIONAL COORDINATES
        CALL READPDB (ONLYCA,OKAY,IPDB,JCRY,R,TM,TMI,XOR,NOR)        
        IF (OKAY )  THEN
         GO TO 7785
        ELSE 
         STOP ' PROBLEMS READING  PDB FILE '
         END IF
        
C . . AS BEFORE BUT FOR A CAMBRIDGE DATA BANK TYPE OF FILE
 7780   OPEN(UNIT=ICAM,NAME=INFIL2,STATUS='OLD',READONLY
     .  ,DEFAULTFILE='MOTIF.CDB')
        WRITE(6,2995)
        READ(5,3000) ANSWER
        IF (ANSWER.EQ. 'Y'.OR. ANSWER.EQ. 'y' ) THEN
          ALLBUTH = .TRUE.
        ELSE
          ALLBUTH = .FALSE.
          WRITE(6,2996)
        END IF
C . .
C . . NOW READ THE ATOMS AND WRITE THEM TO SCRATCH FILE JCRY
C . . IN FRACTIONAL COORDINATES FOR CRYSTAL GENERATION
C . .
        CALL READCDB (ALLBUTH,ICAM,JCRY,R,TM,TMI,XOR,NOR)
C . .
 7785   CONTINUE

C . .
C . . CRYSTAL GENERATION OPTION
C .
        CALL XTAL(JCRY,ICRY,OUTFIL1,OUTFIL4,NEQT)
        ELSE 
C . . END IF NO CRYSTAL GENERATION        
 7020   CONTINUE
        TYPE 7030
 7030   FORMAT(5X, 'NO CRYSTAL GENERATION: I AM DONE! ',/)
        END IF       
        STOP '!VIVE LA REPUBLIQUE!'
        END
C . .

C==============================================================================
      SUBROUTINE ATOMSIN(ALLBUTH,IN,JCRY,NUMATO,ROT,TM,TMI,XOR,NOR)
C . .
C . . ROUTINE TO READ ATOMS MOTIF.MSF FILE
C . . THIS FORMAT IS INTERNAL FOR THE CMD SYSTEM AT ABBOTT LABORATORIES.
C . . 
C . . READS ATOMS APPLIES THE ROTATION AND TRANSLATION
C . . OPERATIONS AND WRITES THEM OUT ON SCRATCH FILE SPG1.CRY
C . . ON LOGICAL UNIT JCRY AS FRACTIONAL COORDINATES FOR THE
C . . SPECIFIED UNIT CELL.
C . .
C . .
C . . THE MATRICES ROT,TM,TMI AND XOR ARE TRANSFERRED FROM THE MAIN PROGRAM:
       DIMENSION ROT(3,4),TM(3,3),TMI(3,3),XOR(3,4)
C . . ARRAYS VI,VO ARE WORKING ARRAYS FOR INPUT AND OUTPUT OF COORDINATES:
       DIMENSION VI(3),VO(3),RF(3),RO(3)
C . .
C . .
       LOGICAL ALLBUTH
C . .
       CHARACTER*6 IATN
       CHARACTER*2 CA
       CHARACTER*1 IATO,H,N,C,O,S
C . .
       DATA H,N,C,O,S/'H','N','C','O','S'/
       DATA CA /'CA'/
       DATA FIN /'END '/
C . .
C . . . .ZERO OUT A FEW VARIABLES
        NUIN = 0
        NUAU = 0
        NID = 0
        NH = 0
        NCA = 0
C . . KEEP TRACK OF OXYGENS IN CASE THERE ARE TOO MANY ATOMS
        NOX=0
C . .EVERY NMOD(THG) ATOM WILL BE PRINTED
        NMOD=10  


        DO 60 IK=1,NUMATO
        READ(IN,8010,end=1000) IATN,X,Y,Z,CHAR,IATO
 8010   FORMAT(A6,1X,4(F10.5,1X),A2,4X)
        NUIN = NUIN+ 1
        IF (IATN(1:2).EQ.CA) NCA= NCA + 1 
        IF ( ALLBUTH) THEN
          IF(IATO.EQ.H)  GO TO 55
        ELSE 
          IF(IATO.EQ.H)  GO TO 55 
          IF(IATO.EQ.O)  GO TO 50
        END IF
C . .          IF (IATN(1:2).NE.CA) GO TO 60
        IF(IATO.EQ.C) ITYPE=1
        IF(IATO.EQ.N) ITYPE=2
        IF(IATO.EQ.O) ITYPE=3
        IF(IATO.EQ.S) ITYPE=4

C .. APPLY ROTATION AND TRANSLATION IN ORTHOGONAL SET OF AXIS
C .. VECTOR POST-MULTIPLIES MATRIX   XX= ROT(I,J)*X
        VI(1) = X
        VI(2) = Y
        VI(3) = Z
        CALL MATVEC1(VO,ROT,VI)
C . . . 
C . . . CHANGE COORDINATES TO FRACTIONAL COORDINATES IN
C ..... OBLIQUE SYSTEM . . . . . . . . .
C . . . . MATRIX TMI (ABOVE) PRE-MULTIPLIES VECTOR 
C . . . . ALGEBRA FROM  NOTES ON MOLECULAR
C . . . . REPLACEMENT FROM J.E. JOHNSON AT PURDUE UNIVERS.
C . . . .  x(fractional) = TMI*X (angstroms orthogonal)
        DO M=1,3
          VI(M) = VO(M)
          RF(M) =0.
          RO(M) = 0.
        END DO
          CALL MATVEC2(RF,TMI,VI)
C . . .
C . . .         
C . . . RF(3) CONTAINS NOW THE FRACTIONAL COORDINATES
C . . . LOOP TO TAKE INTO ACCOUNT THE DIFFERENT BRAVAIS LATTICES
        DO IOR=1,NOR
          NUAU= NUAU + 1
          NID = NUAU +99 
             DO K=1,3
               RF(K) = XOR(K,IOR) + RF(K)
             END DO
C . . WRITE OUT THE TEMPORARY FILE SPG1.CRY CONTAINING
C . . FRACTIONAL COORDINATES OF ATOMS  
          WRITE(JCRY,800)IATN,NID,RF
          IF(MOD(NUIN,NMOD).EQ.0)
     .    WRITE(6,25) X,Y,Z,VO,RF,ITYPE,IATN
 25       FORMAT(2X,9F10.5,I5,1X,A12)
 800      FORMAT(A1,I3,6X,3F10.5)
        END DO 
C ...
        GO TO 60
 50     CONTINUE
        NOX = NOX + 1
        GO TO 60
 55     NH = NH + 1
 60     CONTINUE
 1000   CONTINUE
        IF (NCA.NE.0 ) WRITE(6,6980) NCA
 6980   FORMAT(5X, ' INPUT MSF FILE WAS A PROTEIN FILE WITH ',
     .2X, I5,2X, ' CA POSITIONS ',/)
        WRITE(6,30) NUIN,NUAU,NH,NOX
 30     FORMAT(5X, 'NUMBER OF ATOMS READ IN',2X,I10,/
     .,5X, 'NUMBER OF ATOMS ON ASYM UNIT FILE 2',2X,I10,/
     .,5X, 'NUMBER OF HYDROGEN ATOMS EXCLUDED ',2X,I10,/
     .,5X, 'NUMBER OF OXYGEN ATOMS EXCLUDED   ',2X,I10)
        WRITE(JCRY,6990) FIN
 6990   FORMAT(A)
        CLOSE (UNIT=JCRY)
        RETURN
        END
C . .

C==============================================================================
      SUBROUTINE READCDB(ALLBUTH,IN,JCRY,ROT,TM,TMI,XOR,NOR)
C . .
C . . ROUTINE TO READ ATOMS IN THE SO-CALLED FREE-FORMAT DATA ENTRY
C . . OF THE CAMBRIDGE DATA BANK. (PAG. 11.22) OF THE
C . . CAMBRIDGE CRYSTALLOGRAPHIC DATABASE USER MANUAL (April 1978) 
C . . Such an entry is composed of the following records:
C . . 
C . .    TITLE Record
C . .    CELL  record
C . .    SYMM  record(s)
C . .    Atomic coordinate records
C . .    END  record 
C . .
C . . Since the Cell constants and symmetry operators are already
C . . known from the space group part of the program, only the
C . . atomic records are used in this routine. Coordinates are assumed
C . . to be orthogonal coordinates of a molecular motif placed in some
C . . asymmetric unit of the crystal.
C . .           ============> THEREFORE, <===========
C . . ===> NEITHER CELL, NOR SYMM RECORDS ARE EXPECTED IN THIS FILE <=======
C . .
C . . ROUTINE READS ATOMS, APPLIES THE ROTATION AND TRANSLATION
C . . OPERATIONS AND WRITES THEM OUT ON SCRATCH FILE SPG1.CRY
C . . ON LOGICAL UNIT JCRY AS FRACTIONAL COORDINATES FOR THE
C . . SPECIFIED UNIT CELL.
C . .
C . .
C . . THE MATRICES ROT,TM,TMI AND XOR ARE TRANSFERRED FROM THE MAIN PROGRAM:
       DIMENSION ROT(3,4),TM(3,3),TMI(3,3),XOR(3,4)
C . . ARRAYS VI,VO ARE WORKING ARRAYS FOR INPUT AND OUTPUT OF COORDINATES:
       DIMENSION VI(3),VO(3),RF(3),RO(3)
C . .
C . .
       LOGICAL ALLBUTH
C . .
       CHARACTER*80 TITLE
       CHARACTER*5 ATNAME
       CHARACTER*4 OPTION
       CHARACTER*2 ATOMID
       CHARACTER*1 IATO,H,N,C,O,S
C . .
       DATA H,N,C,O,S/'H','N','C','O','S'/
       DATA FIN /'END '/
C . .
C . . . .ZERO OUT A FEW VARIABLES
        NUIN = 0
        NUAU = 0
        NID = 0
        NH = 0
C . . KEEP TRACK OF OXYGENS IN CASE THERE ARE TOO MANY ATOMS
        NOX=0
C . .EVERY NMOD(THG) ATOM WILL BE PRINTED
        NMOD=10  
        READ(IN,106) TITLE
 106    FORMAT(A)
        IF (TITLE(1:5).EQ.'TITLE') THEN
         WRITE(6,107) TITLE
 107    FORMAT(5X, ' TITLE ON CAMBRIDGE DATA BANK FILE: ',/
     .  5X, A)
        ELSE
         WRITE (6,109) 
 109     FORMAT(5X, 'NO TITLE RECORD FOUND ON C.D.B FILE')
         STOP 'PROBLEMS READING CAMBRIGE DATA BANK FILE'        
        END IF  
 100    READ(IN,1105,END=1110) ATNAME,X,Y,Z
 1105   FORMAT(A,5X,3F10.5)
        IF (ATNAME(1:3).EQ.'END') GO TO 1110
        NUIN = NUIN + 1
        IF (ALLBUTH) THEN
C . . EXCLUDE ONLY HYDROGENS
         IF (ATNAME(1:1).EQ. H ) GO TO 55
        ELSE
C . . EXCLUDE HYDROGENS AND OXYGENS
         IF (ATNAME(1:1).EQ. O ) GO TO 50
         IF (ATNAME(1:1).EQ. H ) GO TO 55
        END IF
C .. APPLY ROTATION AND TRANSLATION IN ORTHOGONAL SET OF AXIS
C .. VECTOR POST-MULTIPLIES MATRIX   XX= ROT(I,J)*X
        VI(1) = X
        VI(2) = Y
        VI(3) = Z
        CALL MATVEC1(VO,ROT,VI)
C . . . 
C . . . CHANGE COORDINATES TO FRACTIONAL COORDINATES IN
C ..... OBLIQUE SYSTEM . . . . . . . . .
C . . . . MATRIX TMI (ABOVE) PRE-MULTIPLIES VECTOR 
C . . . . ALGEBRA FROM  NOTES ON MOLECULAR
C . . . . REPLACEMENT FROM J.E. JOHNSON AT PURDUE UNIVERS.
C . . . .  x(fractional) = TMI*X (angstroms orthogonal)
        DO M=1,3
          VI(M) = VO(M)
          RF(M) =0.
          RO(M) = 0.
        END DO
          CALL MATVEC2(RF,TMI,VI)
C . . .
C . . .         
C . . . RF(3) CONTAINS NOW THE FRACTIONAL COORDINATES
C . . . LOOP TO TAKE INTO ACCOUNT THE DIFFERENT BRAVAIS LATTICES
        DO IOR=1,NOR
          NUAU= NUAU + 1
          NID = NUAU +99 
             DO K=1,3
               RF(K) = XOR(K,IOR) + RF(K)
             END DO
C . . WRITE OUT THE TEMPORARY FILE SPG1.CRY CONTAINING
C . . FRACTIONAL COORDINATES OF ATOMS IN ASSUMMED UNIT CELL.  
          WRITE(JCRY,800)ATNAME,NID,RF
          IF(MOD(NUIN,NMOD).EQ.0)
     .    WRITE(6,25) X,Y,Z,VO,RF,ATNAME
 25       FORMAT(2X,9F10.5,1X,A5)
 800      FORMAT(A1,I3,6X,3F10.5)
        END DO 
C ...
 50     CONTINUE
        NOX = NOX + 1
        GO TO 60
 55     NH = NH + 1
 60     CONTINUE
        GO TO 100
 1110   CONTINUE
        WRITE(6,30) NUIN,NH,NOX
 30     FORMAT(5X, 'NUMBER OF ATOMS READ IN',2X,I10,/
     .,5X, 'NUMBER OF HYDROGEN ATOMS EXCLUDED ',2X,I10,/
     .,5X, 'NUMBER OF OXYGEN ATOMS EXCLUDED   ',2X,I10)
        WRITE(JCRY,6990) FIN
 6990   FORMAT(A)
        CLOSE (UNIT=JCRY)
        RETURN
        END


C==============================================================================
        REAL FUNCTION VOLUME(CELL)
C . . . GIVEN REAL CELL COMPUTES VOLUME 
        DIMENSION CELL(6)
C . . . CELL(4)-CELL(6) CONTAIN COSA,COSB,COSG RESPECTIVELY
C . . . .
        COSA= CELL(4)
        COSB= CELL(5)
        COSG= CELL(6)
        A1 = CELL(1)
        A2 = CELL(2)
        A3 = CELL(3)
        TEMP =(1.-COSA*COSA-COSB*COSB-COSG*COSG+2.*COSA*COSB*COSG)
        VOLUME =A1*A2*A3*SQRT(TEMP)
C . . . .
        RETURN
              END
       

C==============================================================================
      FUNCTION VECSYM(S,TEXT,ICOL) 
C 
C A FORTRAN FUNCTION SUBPROGRAM TO SCAN A SYMMETRY CARD AND BUILD 
C THE CORRESPONDING SYMMETRY VECTOR.  THE CONTENTS OF THE 
C CARD AND ERRORS ENCOUNTERED ARE WRITTEN OFF-LINE IN THE CALLING 
C PROGRAM.  THE FUNCTIONAL VALUE IS NORMALLY ZERO, WHEN AN ERROR HAS
C BEEN DETECTED THE NUMBERS 1-4 ARE GIVEN BACK: 
C ERROR-NUMBER 1:  INVALID CHARACTER AT POSITION .... 
C              2:  SYNTAX ERROR AT NONBLANK POSITION .... 
C              3:  BAD COMMAS 
C              4:  DETERMINANT OF ROTATION MATRIX NOT +1 OR -1
C 
C EXPLANATION OF ARGUMENTS
C S IS THE 1ST ELEMENT OF THE 3 ARRAY WHERE THE VECSYM  IS TO BE 
C STORED. 
C TEXT  IS THE ARRAY OF THE TEXTRECORD WHICH IS TO BE SCANNED  (LENGTH:
C 36 CHARACTERS, CONTIGUOUSLY)
C ICOL POINTS - IN CASE OF ERRORS - TO THE POSITION OF AN INVALID 
C CHARACTER OR AN SYNTAX ERROR WITHIN THE SYMMETRY CARD
C
C *****CHECKED BY C.A-Z FEB.13.87 TO ALLOW FRACTIONS WITH FIVE (5)
C *****IN THE NUMERATOR WHICH OCCUR IN THE POSITION OF SOME OF THE
C *****SYMMETRY OPERATORS IN SPACE GROUP R32 (HEXAGONAL CELL).
C
C ***** UPDATED BY C.A-Z MAY 9.1986 TO BE ABLE TO HANDLE FRACTIONS
C ***** LIKE 1/12 ETC. WHICH OCCUR IN HEXAGONAL SPACE GROUPS.
C ***** INPUT FOR 1/12 SHOULD BE 1/9 TO KEEP ALL DENOMINATORS TO
C ***** A SINGLE CHARACTER. SYMVEC PROGRAM WOULD APPLY THE APPROPIATE
C ***** TRANSLATIONS.
C
C ***** UPDATED BY C.A-Z MAY 6.1986 TO BE ABLE TO HANDLE FRACTIONS
C ***** LIKE 1/8 ETCC. WHICH ARE PRESENT IN THE TETRAGONAL SPACE GROUPS 
C *****
C ***** WRITTEN BY
C *****    C. ABAD-ZAPATERO FROM ORIGINAL MATSYM FUNCTION
C *****    TO GIVE THE VECTOR COORDINATES AND NOT THE SYMMETRY OPERATOR
C *****    **ABBOTT LABORATORIES * DECEMBER 3. 1985.** 
C *****    MODIFIED FOR PDP-11 06/JAN/78 S.J. REMINGTON
C *****     LAST UPDATE 07/10/76   MRS                WRS:PR2.MATSYM
C *****     PREVIOUS UPDATES                       07/07/75  06/11/72
C 
C 
      DIMENSION S(3),ICHAR(36)
      LOGICAL*1 TEXT(36),CHAR,VALID(17),BLANK
      EQUIVALENCE (BLANK,VALID(17))
C . .  THE CHARACTER '9' IN THE DENOMINATOR PLAYS THE ROLE OF '12'
      DATA VALID/'0','1','2','3','4','5','6','8','9',
     .'X','Y','Z','-','+','/',',',' '/
C . .
C 
C----INITIALIZE SIGNALS 
      ICOL = 0
      IFIELD = 0
C----CLEAR SYMOP ARRAY
      DO  I=1,3
        S(I) = 0.0
      END DO
C 
C----TEST THE SYMMETRY TEXT (TEXT) FOR ILLEGAL CHARACTERS
C 
C      THE STATEMENT IS SCANNED AND ARRAY ICHAR IS FILLED WITH THE INDICES
C      OF EACH ELEMENT OF TEXT INTO THE ARRAY 'VALID' (I.E. IF THE SEVENTH
C      CHAR OF 'TEXT' -IGNORING BLANKS- IS 'Y', ICHAR(7)=8)
C      VALID CHARACTERS ARE, FOR VECSYM:012345689XYZ-+/,   .(17)
C . . .THE ONLY CHANGE NECESSARY FOR CRYSTALLOGRAPHIC FRACTIONS
C . . .IS TO REPLACE 12, IN THE DENOMINATOR BY 9. THE CODE MAKES IT 1/12
C . . .BUT IN THIS FORM ALL FRACTIONS REQUIRE ONLY ON CHARACTER FOR THE
C . . .NUMERATOR AND ANOTHER FOR THE DENOMINATOR.
C . . .
C . . .THE CHARACTER FIVE IS FORBIDEN IN THE DENOMINATOR FOR
C . . .THE TRANSLATIONAL COMPONENT  
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      ICOLMX = 0
      DO 20 I=1,36
      CHAR = TEXT(I)
      IF(CHAR .EQ. BLANK) GO TO 20
      ICOLMX = ICOLMX + 1
      DO 10 J=1,16
      IF( VALID(J) .NE. CHAR) GO TO 10
      ICHAR(ICOLMX) = J
      GO TO 20
   10 CONTINUE
C--IF GOT TO HERE, THE CHAR ISN'T IN 'VALID'
      JTXTSC = I
      GO TO 9975
   20 CONTINUE
C . . DEBUG 
C      WRITE(6,3000) ICHAR
C3000 FORMAT(36(1X,I2))
C
C----BEGIN FIELD LOOP 
C 
 101  IFIELD = IFIELD + 1 
      IF(IFIELD-3) 104,104,103
C  HAVE ALL CHARACTERS  BEEN ANALYSED IN 3 FIELDS?
 103  IF(ICOL-ICOLMX) 9987,1000,1000
C---NO MORE THAN THREE FIELDS PERMITTED 
 104  SIGN = 1.0
      ICOL = ICOL+1 
C----MARCH ON, FIND WHATS NEXT
      IFNUM=ICHAR(ICOL) 
      IF (14-IFNUM) 9980,110,110
C----TEST FOR SIGNS 
 110  IF (13-IFNUM) 118,115,112 
C----TEST TO DISTINGUISH BETWEEN NUMBERS AND LETTERS
 112  IF (9-IFNUM) 125,130,130
C----FOUND A SIGN 
C     COME HERE FOR MINUS 
 115  SIGN = -1.0 
C     COME HERE FOR PLUS
 118  ICOL = ICOL + 1 
C----NEXT CHARACTER MUST BE A NUMBER OR LETTER
      IFNUM=ICHAR(ICOL) 
      IF (IFNUM-12) 122,122,9980 
 122  IF (IFNUM-9) 130,130,125
C----A LETTER 
 125  IFNUM = IFNUM - 9 
      S(IFIELD) = SIGN
      GO TO 150 
C----A NUMBER, FLOAT IT INTO DIGIT
C----NUMERATOR OF FRACTION
 130  DIGIT = IFNUM - 1
C----
C----
      ICOL = ICOL + 1 
C----IS NEXT CHARACTER A SLASH
      IF(ICHAR(ICOL).EQ.15)  GO TO 135
C----NO SLASH, TRANSLATION TERM COMPLETE
      ICOL=ICOL-1 
      GO TO  140
C----A SLASH IS NEXT, IS THERE A NUMBER AFTER IT
 135  ICOL = ICOL + 1 
      IFNUM=ICHAR(ICOL) 
C----THE FLOATING IS DIFFERENT DEPENDING WHETHER 
C----NUMBER IS 0-6 , 8, 9 OR GREATER 
      IF (IFNUM-8) 137,138,139 
C----MAKE FRACTION;   '5' NOT ALLOWED IN DENOMINATOR
C----NUMBERS 0,6  DIGIT = IFNUM - 1
C----
  137 IF(IFNUM.EQ.6) GO TO 9980    
      DIGIT = DIGIT/FLOAT(IFNUM-1)
      GO TO 140
C----NUMBER 8  SPECIAL CASE
C----  8 ===>  DIGIT =  DIGIT/FLOAT(IFNUM) 
  138 DIGIT = DIGIT/FLOAT(IFNUM)   
      GO TO 140
C----NUMBER 9 SPECIAL CASE
C----  8 ===>  DIGIT = IFNUM + 3
  139 IF (IFNUM.GT.9 ) GO TO 9980
      DIGIT = DIGIT/FLOAT(IFNUM+3)
C---- 
C----ACCUMULATE THE VECTOR COMPONENT
 140  S(IFIELD) = S(IFIELD) + SIGN*DIGIT
C----TEST TO SEE IF THERE ARE MORE CHARACTERS IN THIS SUBFIELD
C    OR A NEW SUBFIELD
 150  ICOL = ICOL + 1 
      SIGN = 1.0
C----TEST IF ALL DONE 
      IF(ICOL - ICOLMX) 151,151,1000
 151  CONTINUE
C----A SUBFIELD BEGINS WITH A PLUS OR MINUS 
      IF(ICHAR(ICOL).EQ.14)  GO TO 118
      IF(ICHAR(ICOL).EQ.13)  GO TO 115
C----A COMMA MUST BE NEXT UNLESS ICOL IS ICOLMX+1 WHICH MEANS THE END 
      IF(ICHAR(ICOL).EQ.16)   GO TO 101 
 163  IF (ICOLMX+1-ICOL) 9980,1000,9980 
C----EVERYTHING SEEMS OK  
 1000 CONTINUE
      IF(IFIELD - 3) 9987,1002,9987 
 1002 CONTINUE
C----LEGAL RETURN 
      VECSYM = 0
 2005 RETURN
C 
C 
C----ILLEGAL CHARACTER ON SYMMETRY CARD 
 9975 VECSYM = 1
      ICOL = JTXTSC 
      GO TO 2005
C 
C----SYNTAX ERROR AT NONBLANK CHARACTER ICOL
 9980 VECSYM = 2
      GO TO 2005
C 
C----INCORRECT NUMBER OF COMMAS 
 9987 VECSYM = 3
      GO TO 2005
      END 


C==============================================================================
        FUNCTION MATSYM(S,TEXT,ICOL) 
C 
C A FORTRAN FUNCTION SUBPROGRAM TO SCAN A SYMMETRY CARD AND BUILD 
C THE CORRESPONDING SYMMETRY-OPERATION MATRIX.  THE CONTENTS OF THE 
C CARD AND ERRORS ENCOUNTERED ARE WRITTEN OFF-LINE IN THE CALLING 
C PROGRAM.  THE FUNCTIONAL VALUE IS NORMALLY ZERO, WHEN AN ERROR HAS
C BEEN DETECTED THE NUMBERS 1-4 ARE GIVEN BACK: 
C ERROR-NUMBER 1:  INVALID CHARACTER AT POSITION .... 
C              2:  SYNTAX ERROR AT NONBLANK POSITION .... 
C              3:  BAD COMMAS 
C              4:  DETERMINANT OF ROTATION MATRIX NOT +1 OR -1
C 
C EXPLANATION OF ARGUMENTS
C S IS THE 1ST ELEMENT OF THE 3X4 ARRAY WHERE THE SYM. OP. IS TO BE 
C STORED. 
C TEXT  IS THE ARRAY OF THE TEXTRECORD WHICH IS TO BE SCANNED  (LENGTH:
C 36 CHARACTERS, CONTIGUOUSLY)
C ICOL POINTS - IN CASE OF ERRORS - TO THE POSITION OF AN INVALID 
C CHARACTER OR AN SYNTAX ERROR WITHIN THE SYMMETRY CARD 
C
C ***** MODIFIED FOR TABLES PROGRAM BY C.ABAD-ZAPATERO
C ***** ABBOTT LABORATORIES . DECEMBER 6. 1985. 
C 
C     MODIFIED FOR PDP-11 06/JAN/78 S.J. REMINGTON
C ***** LAST UPDATE 07/10/76   MRS                WRS:PR2.MATSYM
C ***** PREVIOUS UPDATES                              07/07/75  06/11/72
C 
C 
      DIMENSION S(3,4), ICHAR(36)
      LOGICAL*1 TEXT(36),CHAR,VALID(14),BLANK
      EQUIVALENCE (BLANK,VALID(14))
      DATA VALID/'1','2','3','4','5','6','X','Y','Z','-','+','/',',',
     . ' '/
C 
C----INITIALIZE SIGNALS 
      ICOL = 0
      IFIELD = 0
C----CLEAR SYMOP ARRAY
      DO 4 J=1,4
      DO 4 I=1,3
 4    S(I,J) = 0.0
C 
C----TEST THE SYMMETRY TEXT (TEXT) FOR ILLEGAL CHARACTERS
C 
C      THE STATEMENT IS SCANNED AND ARRAY ICHAR IS FILLED WITH THE INDICES
C      OF EACH ELEMENT OF TEXT INTO THE ARRAY 'VALID' (I.E. IF THE SEVENTH
C      CHAR OF 'TEXT' -IGNORING BLANKS- IS 'Y', ICHAR(7)=8)
C
      ICOLMX = 0
      DO 20 I=1,36
      CHAR = TEXT(I)
      IF(CHAR .EQ. BLANK) GO TO 20
      ICOLMX = ICOLMX + 1
      DO 10 J=1,13
      IF( VALID(J) .NE. CHAR) GO TO 10
      ICHAR(ICOLMX) = J
      GO TO 20
   10 CONTINUE
C--IF GOT TO HERE, THE CHAR ISN'T IN 'VALID'
      JTXTSC = I
      GO TO 9975
   20 CONTINUE
C
C----BEGIN FIELD LOOP 
C 
 101  IFIELD = IFIELD + 1 
      IF(IFIELD-3) 104,104,103
C  HAVE ALL CHARACTERS  BEEN ANALYSED IN 3 FIELDS?
 103  IF(ICOL-ICOLMX) 9987,1000,1000
C---NO MORE THAN THREE FIELDS PERMITTED 
 104  SIGN = 1.0
      ICOL = ICOL+1 
C----MARCH ON, FIND WHATS NEXT
      IFNUM=ICHAR(ICOL) 
      IF (11-IFNUM) 9980,110,110
C----TEST FOR SIGNS 
 110  IF (10-IFNUM) 118,115,112 
C----TEST TO DISTINGUTSH BETWEEN NUMBERS AND LETTERS
 112  IF (6-IFNUM) 125,130,130
C----FOUND A SIGN 
C     COME HERE FOR MINUS 
 115  SIGN = -1.0 
C     COME HERE FOR PLUS
 118  ICOL = ICOL + 1 
C----NEXT CHARACTER MUST BE A NUMBER OR LETTER
      IFNUM=ICHAR(ICOL) 
      IF (IFNUM-9) 122,122,9980 
 122  IF (IFNUM-6) 130,130,125
C----A LETTER 
 125  IFNUM = IFNUM - 6 
      S(IFIELD,IFNUM) = SIGN
      GO TO 150 
C----A NUMBER, FLOAT IT INTO DIGIT
 130  DIGIT = IFNUM 
      ICOL = ICOL + 1 
C----IS NEXT CHARACTER A SLASH
      IF(ICHAR(ICOL).EQ.12)  GO TO 135
C----NO SLASH, TRANSLATION TERM COMPLETE
      ICOL=ICOL-1 
      GO TO  140
C----A SLASH IS NEXT, IS THERE A NUMBER AFTER IT
 135  ICOL = ICOL + 1 
      IFNUM=ICHAR(ICOL) 
      IF (IFNUM-6) 138,138,9980 
C----MAKE FRACTION;   '5' NOT ALLOWED IN DENOMINATOR
  138 IF(IFNUM.EQ.5) GO TO 9980 
      DIGIT = DIGIT/FLOAT(IFNUM)
C----ACCUMULATE THE VECTOR COMPONENT
 140  S(IFIELD,4) = S(IFIELD,4) + SIGN*DIGIT
C----TEST TO SEE IF THERE ARE MORE CHARACTERS IN THIS SUBFIELD
C    OR A NEW SUBFIELD
 150  ICOL = ICOL + 1 
      SIGN = 1.0
C----TEST IF ALL DONE 
      IF(ICOL - ICOLMX) 151,151,1000
 151  CONTINUE
C----A SUBFIELD BEGINS WITH A PLUS OR MINUS 
      IF(ICHAR(ICOL).EQ.11)  GO TO 118
      IF(ICHAR(ICOL).EQ.10)  GO TO 115
C----A COMMA MUST BE NEXT UNLESS ICOL IS ICOLMX+1 WHICH MEANS THE END 
      IF(ICHAR(ICOL).EQ.13)   GO TO 101 
 163  IF (ICOLMX+1-ICOL) 9980,1000,9980 
C----EVERYTHING SEEMS OK SEE IF DETERMINATE IS A NICE + OR - 1. 
 1000 CONTINUE
      D=S(1,1)*S(2,2)*S(3,3)
     $ -S(1,2)*S(2,1)*S(3,3)
     $ -S(1,1)*S(2,3)*S(3,2)
     $ -S(1,3)*S(2,2)*S(3,1)
     $ +S(1,2)*S(2,3)*S(3,1)
     $ +S(1,3)*S(2,1)*S(3,2)
      IF(ABS(ABS(D) - 1.0) -.0001) 1001,9985,9985 
 1001 CONTINUE
      IF(IFIELD - 3) 9987,1002,9987 
 1002 CONTINUE
C----LEGAL RETURN 
      MATSYM = 0
 2005 RETURN
C 
C 
C----ILLEGAL CHARACTER ON SYMMETRY CARD 
 9975 MATSYM = 1
      ICOL = JTXTSC 
      GO TO 2005
C 
C----SYNTAX ERROR AT NONBLANK CHARACTER ICOL
 9980 MATSYM = 2
      GO TO 2005
C 
C----DETERMINANT IS NOT + OR - 1. 
 9985 MATSYM = 4
      GO TO 2005
C 
C----INCORRECT NUMBER OF COMMAS 
 9987 MATSYM = 3
      GO TO 2005
      END 

C==============================================================================
        SUBROUTINE TMATRIX(CELL,VOL,T)
C. . . . SETS UP THE ORTHOGONALIZATION MATRIX TO GO FROM
C. . . . FRACTIONAL COORDINATES IN x,y,z (CRYSTAL SPACE)
C. . . . TO AN ORTHOGONAL FRAME OF REFERENCE X,Y,Z 
C. . . . NEW AXIS WILL BE:
C        X' ALONG A-AXIS, Y' IN THE (A,B) PLANE AND Z' ALONG C*.
        DIMENSION CELL(6),T(3,3)
        A = CELL(1)
        B = CELL(2)
        C = CELL(3)
        COSA= CELL(4)
        COSB= CELL(5)
        COSG= CELL(6)
        SINA = SQRT(1.-COSA*COSA)
        SINB = SQRT(1.-COSB*COSB)
        SING = SQRT(1.-COSG*COSG)
        T(1,1) = A
        T(1,2) = B*COSG
        T(1,3) = C*COSB
        T(2,1) = 0.
        T(2,2) = B*SING
        T(2,3)=C*(COSA-COSB*COSG)/SING
        T(3,1) = 0.
        T(3,2) = 0.
        T(3,3) = VOL/(A*B*SING)
        RETURN
        END

C==============================================================================
       SUBROUTINE INVERS(A,AI)
C ...  TO COMPUTE THE INVERS OF MATRIX A: AI
C ...  ONLY VALID FOR UPPER TRIANGULAR MATRICES
       DIMENSION A(3,3),AI(3,3)
C .
c . . . FOR UPPER TRIANGULAR MATRIX DET = a(1,1)*a(2,2)*a(3,3)
                DET = A(1,1)*A(2,2)*A(3,3)
C . . .
		AI(1,1)=1./A(1,1)
                AI(1,2) = -A(1,2)/(A(1,1)*A(2,2))
		AI(1,3) = (A(1,2)*A(2,3)-A(1,3)*A(2,2))/DET
		AI(2,1) = 0.0
		AI(2,2) = 1./A(2,2)
		AI(2,3) = -A(1,1)*A(2,3)/DET
		AI(3,1) = 0.0
		AI(3,2) = 0.0
		AI(3,3) = 1./A(3,3)
		RETURN
		END


C==============================================================================
    		SUBROUTINE ROTMAT(CAPA,PSI,PHI,R)
C . . ROUTINE TO FORM THE ROTATION MATRIX R(I,J) ABOUT
C . . AN AXIS DEFINED BY PSI,PHI (POLAR COORDINATES) AND CAPA
C . . DEGREES
C . . INTL. TABLES VOL II  PAG. 63
		DIMENSION R(3,3)
		REAL*8 L1,L2,L3
		L1 = COS(PHI)*SIN(PSI)
		L2 = SIN(PHI)*SIN(PSI)
		L3 = COS(PSI)
		CA = COS(CAPA)
		SA = SIN(CAPA)
		CAM1 = 1. - CA
		R(1,1) = CA + L1*L1*CAM1
		R(1,2) = L1*L2*CAM1+L3*SA
		R(1,3) = L1*L3*CAM1-L2*SA
		R(2,1) = L2*L1*CAM1-L3*SA
		R(2,2) = CA + L2*L2*CAM1
		R(2,3) = L2*L3*CAM1+L1*SA
		R(3,1) = L3*L1*CAM1+L2*SA
		R(3,2) = L3*L2*CAM1-L1*SA
		R(3,3) = CA + L3*L3*CAM1
		RETURN
		END
C . . 

C==============================================================================
        SUBROUTINE HEADMSF(IO,LIN1,LIN2,DATE,TIME,PROG,VERS,IAB,
     .             NATS,NBONS)
C . .  READS THE HEADER OF A .MSF TYPE FILE ON UNIT IO
C . .  AND RETURNS 
C . .  THE IMPORTANT VALUES: NUMBER OF ATOMS (NA) 
C . .                      : NUMBER OF BONDS (NB)
C . . . THE ***.MSF FILE CONTAINS 3 HEADER RECORDS DESCRIBED AS
C . . . TYPE       FORMAT       FIELD  (FROM FRED THOMAS)
C . .     1          80A1         1    TITLE
C . .    
C . .    2A          80A          1    SUBTITLE RECORD
C . .
C . .
C . .    2B         A8,1X        1    DATE OF CREATION, OBJECT FILE
C . .               A8,1X        2    TIME OF CREATION, OBJECT FILE
C . .               A9,1X        4    CREATING PROGRAM OBJECT FILE
C . .               A9,1X        5    VERSION OF PROGRAM 
C . .                I6,1X        5    ABBOTT  NUMBER
C . .                I5,1X        6    NUMBER OF ATOMS
C . .                I5           7    NUMBER OF BONDS(ORBITALS)
C . .
C . .      FOLLOWING RECORDS
C . .
C . .  6A1,1X,4(F10.5,1X),A6;  ATOM NAME,X,Y,Z,CHARGE,ATOMIC SYMBOL
C . .
        CHARACTER*80 LIN1,LIN2
        CHARACTER*9 PROG,VERS
        CHARACTER*8 DATE,TIME

        READ(IO,10) LIN1
        WRITE(6,15) LIN1
 15     FORMAT(5X,'HEADER TITLE IN INPUT FILE.MSF ',/,5X,A)        
        READ(IO,10) LIN2
        WRITE(6,20) LIN2
 20     FORMAT(5X,'HEADER SUBTITLE IN FILE.MSF ',/,5X,A )
        READ(IO,25) DATE,TIME,PROG,VERS,IAB,NA,NB
 25     FORMAT(A8,1X,A8,1X,A9,1X,A9,1X,I6,1X,I5,1X,I5)
        WRITE(6,30) DATE,TIME,PROGR,VERS,IAB,NA,NB 
 30     FORMAT(5X, 'DATE FILE CREATION:  ', 1X,A8,/
     .      5X, 'TIME OF FILE CREATION:  ', 1X,A8,/
     .      5X, 'CREATING PROGRAM : ',A9,1X,/
     .      5X, 'VERSION OF CREATING PROGRAM ',A9,1X,/
     .      5X, 'ABBOTT COMPOUND NUMBER ', 1X,I6,/
     .      5X, 'NUMBER OF ATOMS IN FILE ',1X,I5,/
     .      5X, 'NUMBER OF BONDS IN FILE ',1X,I5,//)
 10     FORMAT(A)
        NATS = NA
        NBONS= NB
        RETURN
        END
C . .

C==============================================================================
        SUBROUTINE OUTCELL
     .  (IO,LIN1,LIN2,DAT,TIM,PROG,VERS,IAB,NA,NB)
C . .  WRITES THE CELL.MSF FILE ON UNIT IO
C . .  AND RETURNS
C . .  BASED ON TRANSFERRED INFORMATION IN CALLING ROUTINE
C . . 
C . .  THE IMPORTANT VALUES: NUMBER OF VERTICES (ATOMS): NA 
C . .                      : NUMBER OF LINES (BONDS): NB
C . . . THE ***.MSF FILE CONTAINS 3 HEADER RECORDS DESCRIBED 
C . . . IN ROUTINE HEADMSF (PREVIOUS ROUTINE)
C . . .
C . . .  FOLLOWING RECORDS
C . . .  COORDINATE RECORDS WITH FORMAT:
C . .
C . .  6A1,1X,4(F10.5,1X),A6;  ATOM NAME,X,Y,Z,CHARGE,ATOMIC SYMBOL
C . .
      CHARACTER*1 C
      CHARACTER*4 BONDS
      INTEGER*2 IBONDS
C . .
C . . SET COMMON BLOCKS FOR THE CELL FILE:
C . . 
C . . 
      COMMON /CELL/V000(3,8),BONDS(2,12),IBONDS(2,12),C,WG1,BONTY1     
C . .
        CHARACTER*80 LIN1,LIN2
        CHARACTER*9 PROG,VERS
        CHARACTER*8 DAT,TIM

        WRITE(IO,10) LIN1
        WRITE(6,15) LIN1
 15     FORMAT(5X,'HEADER TITLE IN CELL.MSF ',/,5X,A)        
        WRITE(IO,10) LIN2
        WRITE(6,20) LIN2
 20     FORMAT(5X,'HEADER SUBTITLE IN CELL.MSF ',/,5X,A )
        WRITE(IO,25) DATE,TIME,PROG,VERS,IAB,NA,NB
 25     FORMAT(A8,1X,A8,1X,A9,1X,A9,1X,I6,1X,I5,1X,I5)
        WRITE(6,30) DAT,TIM,PROG,VERS,IAB,NA,NB 
 30     FORMAT(5X, 'DATE FILE CREATION:  ', 1X,A8,/
     .      5X, 'TIME OF FILE CREATION:  ', 1X,A8,/
     .      5X, 'CREATING PROGRAM : ',A9,1X,/
     .      5X, 'VERSION OF CREATING PROGRAM ',A9,1X,/
     .      5X, 'ABBOTT COMPOUND NUMBER ', 1X,I6,/
     .      5X, 'NUMBER OF ATOMS IN FILE ',1X,I5,/
     .      5X, 'NUMBER OF BONDS IN FILE ',1X,I5,//)
 10     FORMAT(A)
C . . WRITE NOW THE VERTICES OF THE UNIT CELL AS CARBON ATOMS
C . . IN THE .MSF TYPE OF FORMAT:
        DO K=1,NA
         ID = 99 + K
         WRITE(IO,8000) C,ID,(V000(JJ,K),JJ=1,3),WG,C
 8000    FORMAT (A,I3,3X,4(F10.5,1X),A2,4X)
        END DO
C . . NOW THE LINES TO MAKE THE PARALLELEPIPED AS PSEUDO BONDS:
         DO IB=1,NB
         WRITE(IO,8005) (BONDS(K,IB), K=1,2), BONTYP
 8005    FORMAT(2(A4,2X),1X,F3.1)
         END DO
C . . ALL DONE! 
        RETURN
        END
C . .

C==============================================================================
        SUBROUTINE OUTOPS
     .  (IO,LIN1,LIN2,DAT,TIM,PROG,VERS,IAB,NSYP,NSYV)
C . .  WRITES THE THE SYMOPS.MSF TYPE FILE ON UNIT IO:
C . .  BASED ON TRANSFERRED INFORMATION IN CALLING ROUTINE
C . .  AND COMMON BLOCK SYMOPS.
C . .  THE IMPORTANT VALUES: NUMBER OF SYMMETRY POINTS (ATOMS): NSYP 
C . .                      : NUMBER OF SYMMETRY VECTORS NSYV(BONDS):NSYV
C . . . THE ***.MSF FILE CONTAINS 3 HEADER RECORDS DESCRIBED 
C . . . IN ROUTINE HEADMSF (PREVIOUS ROUTINE)
C . . .
C . . .  FOLLOWING RECORDS
C . . .  COORDINATE RECORDS WITH FORMAT:
C . .
C . .  6A1,1X,4(F10.5,1X),A6;  ATOM NAME,X,Y,Z,CHARGE,ATOMIC SYMBOL
C . .
      PARAMETER MAXATOM=2000
      PARAMETER MAXSYM =32
      PARAMETER MAXOPS =200
C . .
        CHARACTER*1 IATID
C . .
C . . SET COMMON BLOCKS FOR THE SYMOPS FILE : MAXOPS = 200
C . . 
      COMMON /SYMOPS/SVEC(3,MAXOPS),IATID(MAXOPS),WG2,BONTY2
C . .
        CHARACTER*80 LIN1,LIN2
        CHARACTER*9 PROG,VERS
        CHARACTER*8 DAT,TIM

        WRITE(IO,10) LIN1
        WRITE(6,15) LIN1
 15     FORMAT(5X,'HEADER TITLE IN SYMOPS.MSF ',/,5X,A)        
        WRITE(IO,10) LIN2
        WRITE(6,20) LIN2
 20     FORMAT(5X,'HEADER SUBTITLE IN SYMOPS.MSF ',/,5X,A )
        WRITE(IO,25) DATE,TIME,PROG,VERS,IAB,NSYP,NSYV
 25     FORMAT(A8,1X,A8,1X,A9,1X,A9,1X,I6,1X,I5,1X,I5)
        WRITE(6,30) DAT,TIM,PROG,VERS,IAB,NSYP,NSYV 
 30     FORMAT(5X, 'DATE FILE CREATION:  ', 1X,A8,/
     .      5X, 'TIME OF FILE CREATION:  ', 1X,A8,/
     .      5X, 'CREATING PROGRAM : ',A9,1X,/
     .      5X, 'VERSION OF CREATING PROGRAM ',A9,1X,/
     .      5X, 'ABBOTT COMPOUND NUMBER ', 1X,I6,/
     .      5X, 'NUMBER OF ATOMS IN FILE ',1X,I5,/
     .      5X, 'NUMBER OF BONDS IN FILE ',1X,I5,//)
 10     FORMAT(A)
C . .
C . . . WRITE THE PSEUDO-ATOMS DEFINING THE POSITIONS OF THE
C . . . SYMMETRY ELEMENTS  
        DO ISP=1,NSYP
        J = ISP
        ID = 99+ISP
        WRITE(IO,8000)IATID(J),ID,(SVEC(K,ISP),K=1,3),WG2,IATID(J)
 8000   FORMAT (A,I3,3X,4(F10.5,1X),A,4X)
        END DO
C . . . WRITE THE PSEUDO-BONDS DEFINING THE DIRECTIONS OF THE
C . . . SYMMETRY ELEMENTS IN THE CELL
          DO NV=1,NSYV
           I1 = 1 + (NV-1)*2
           I2 = I1 + 1
           N1 = 100 + (NV-1)*2
           N2 = N1 + 1 
            WRITE(IO,8215) IATID(I1),N1,IATID(I2),N2,BONTY2
 8215       FORMAT(A,I3,2X,A,I3,3X,F3.1)    
          END DO
C . . ALL DONE! 
        RETURN
        END
C . .

C==============================================================================
	SUBROUTINE XTAL(IUNIT,OUNIT,IFILE,OFILE,NSYMM)
	INTEGER IUNIT,OUNIT,NSYMM
	CHARACTER*64 OFILE,IFILE
C
C	SUBROUTINE TO GENERATE THE CRYSTAL STRUCTURE FROM THE
C		UNIT CELL PARAMS AND THE SYMMETRY OPERATORS.
C	REQUIRES THAT IUNIT IS ALREADY OPEN, WILL OPEN OUNIT.
C
C . . THIS ROUTINE RE-COMPUTES SOME GEOMETRICAL PARAMETERS OF THE
C . . CRYSTAL GEOMETRY (FOR INSTANCE VOLUME, AND ORTHO-DEORTHO MATRICES)
C . .
C . . ROUTINES CALLED BY THIS SUBROUTINE ARE PREFIXED BY A TJ (TJMATRIX).
C . . ONE COULD CLEAN THE PROGRAM UP, BUT SINCE IT IS WORKING NOW
C . . I (C. A-Z) DECIDED NOT TO CHANGE IT. MOREOVER, THE ROUTINE XTAL,
C . . WITH THE SUBROUTINES THAT IT CALLS, COULD BE MADE A STAND-ALONE
C . . CRYSTAL GENERATING PROGRAM.  THIS IS THE WAY THAT T.J. ODONNELL
C . . WROTE IT ORIGINALLY. LATER WAS INCORPORATED INTO THE T A B L E S
C . . PROGRAM AS AN OPTION.
C . .
	PARAMETER MAXATOM=2000
	PARAMETER MAXSYM=32
	DIMENSION TMAT(3,3),TTR(3),UTR(3),UNIT(3,3)
	DIMENSION SYMM(3,3,MAXSYM),STR(3,MAXSYM)
	DIMENSION X(3,MAXATOM),UX(3,MAXATOM*MAXSYM)
	DIMENSION AX(3),BX(3),CX(3,MAXATOM)
	DIMENSION UCORN(3),CCORN(3)
	REAL NBC(7),CELL(6)
	EXTERNAL VOL
	CHARACTER*80 TITLE
	CHARACTER*120 FULLTITLE
	CHARACTER*4 OPTION
	CHARACTER*5 ATNAME(MAXATOM)
	LOGICAL ALL27,DCHECK,BOUNDS,WRITEM,SKIP,CWRIT,HEAVY,SKIPATOM
	CHARACTER*1 ANSWER
	CHARACTER*2 ATOMID
	INTEGER TRIM
	REAL XMAX,YMAX,ZMAX,XMIN,YMIN,AMIN,XCT,YCT,ZCT
	EXTERNAL ATOMID,TRIM
	DATA PI/3.141593/
	DATA UNIT /1.,0.,0.,  0.,1.,0.,   0.,0.,1./
	DATA NBC /0.,1.,-1.,2.,-2.,-3.,3./
C
	OPEN(UNIT=IUNIT,NAME=IFILE,STATUS='OLD'
     .,DEFAULTFILE='SPG1.CRY')
	OPEN(UNIT=OUNIT,NAME=OFILE,STATUS='NEW'
     .,DEFAULTFILE='SPG2.CRY')
	READ (IUNIT,106)TITLE
 106	FORMAT(A)
	IF(TITLE(1:5).EQ.'TITLE')THEN
		FULLTITLE=TITLE(6:)
	ELSE
		FULLTITLE=TITLE
	END IF
	WRITE(6,*)FULLTITLE
	WRITE(6,111)
 111	FORMAT(/1X,'YOU HAVE THE FOLLOWING OPTIONS:'/
     &	1X,'1: OUTPUT JUST THE UNIT CELL.'/
     &	1X,'2: OUTPUT A 3x3 ARRAY OF UNIT CELLS.'/
     &	1X,'3: OUTPUT ONLY MOLECULES IN NEIGHBORING CELLS'/
     &	1X,'   THAT ARE WITHIN A GIVEN DISTANCE.'/
     &	1X,'PLEASE ENTER THE OPTION NUMBER YOU''D LIKE: '$)
	READ(5,112)IOPT
	IF(IOPT.EQ.1)THEN
		ALL27=.FALSE.
		DCHECK=.FALSE.
		FULLTITLE=FULLTITLE(:TRIM(FULLTITLE))//' UNIT CELL'
	ELSE IF(IOPT.EQ.2)THEN
		ALL27=.TRUE.
		DCHECK=.FALSE.
		FULLTITLE=FULLTITLE(:TRIM(FULLTITLE))//' 3X3 UNIT CELL NEIGHBORS'
	ELSE IF(IOPT.EQ.3)THEN
		DCHECK=.TRUE.
		ALL27=.TRUE.
		WRITE(6,117)
		READ(5,119)DMAX
		D2MAX=DMAX*DMAX
		WRITE(6,121)
		READ(5,114)ANSWER
		IF(ANSWER.EQ.'Y' .OR. ANSWER.EQ.'y')THEN
			HEAVY=.TRUE.
			FULLTITLE=FULLTITLE(:TRIM(FULLTITLE))//' HEAVY'
		ELSE
			HEAVY=.FALSE.
		END IF
		WRITE(FULLTITLE(TRIM(FULLTITLE)+1:),123)
     &		' ATOMS WITHIN ',DMAX,' ANGSTROM'
 123		FORMAT(A,F5.2,A)
 124		FORMAT(A)
	END IF
	WRITE(6,113)
 113	FORMAT(/1X,'MAY I SHIFT ATOMS WHICH ARE OUTSIDE THE BOUNDARIES '
     &	/1X,'OF THE UNIT CELL SO THAT THEY ARE INSIDE (THIS MAY RESULT'
     &	/1X,'IN FRAGMENTATION OF MOLECULAR ENTITIES)? '$)
	READ(5,114)ANSWER
	IF(ANSWER.EQ.'Y' .OR. ANSWER.EQ.'y')THEN
		BOUNDS=.TRUE.
		FULLTITLE=FULLTITLE(:TRIM(FULLTITLE))//', SHIFTED'
	ELSE
		BOUNDS=.FALSE.
	END IF
 112	FORMAT(I)
 114	FORMAT(A)
 117	FORMAT(/1X,'WHAT IS THE MINIMUM DISTANCE FROM THE UNIT CELL'
     &	/1X,'FOR A MOLECULE TO BE RETAINED? '$)
 119	FORMAT(F)
 121	FORMAT(/1X,'SHOULD THIS DISTANCE APPLY ONLY TO HEAVY ATOMS? '$)
C
	WRITE(6,*)FULLTITLE
	IF(TRIM(FULLTITLE).GT.80)THEN
		WRITE(6,*)'TITLE TRIMMED TO 80 CHARACTERS'
		WRITE(OUNIT,106)FULLTITLE(:80)
	ELSE
		WRITE(OUNIT,106)FULLTITLE(:TRIM(FULLTITLE))
	END IF
C
	PI180=PI/180.
	READ (IUNIT,101) OPTION,A,B,C,ALPHA,BETA,GAMMA
 101	FORMAT(A,4X,6F8.3)
C	IF( ALPHA.NE.90. .OR. GAMMA.NE.90.)STOP 'ONLY MONCLINIC, PLEASE!'
C
C	READ IN THE SYMMETRY ELEMENTS
	DO I=1,NSYMM
	READ (IUNIT,102)OPTION,( ( SYMM(J,K,I),K=1,3 ),STR(J,I),J=1,3   )
C	WRITE (6,3) ((SYMM(J,K,I),K=1,3),J=1,3)
C	WRITE (6,3) (STR(J,I),J=1,3)
	END DO
 3	FORMAT (6F10.5)
 102	FORMAT(A,4X,3(3F4.0,1X,F7.5,1X))
C
C	SET UP THE MATRIX
C
	CELL(1)=A
	CELL(2)=B
	CELL(3)=C
	CELL(4)=COS(PI/180. * ALPHA)
	CELL(5)=COS(PI/180. * BETA )
	CELL(6)=COS(PI/180. * GAMMA)
	CALL TJMATRIX(CELL,VOL(CELL),TMAT)
C
C	  DUMMY, USED FOR ZERO TRANSLATIONS IN MATMPY CALLS
C
	TTR(1)=0.0
	TTR(2)=0.0
	TTR(3)=0.0
C
C	  IF CHECKING DISTANCES, FIND LENGTH OF UNIT CELL DIAGONAL FOR
C	  LATER ESTIMATES OF CLOSENESS OF OTHER TRANSLATED CELLS.
C
	IF(DCHECK)THEN
		UCORN(1)=1.
		UCORN(2)=1.
		UCORN(3)=1.
		CALL MATMPY(TMAT,UCORN,CCORN,TTR)
		DIAG2=2.* (CCORN(1)**2 + CCORN(2)**2 + CCORN(3)**2)
		WRITE(6,24)SQRT(DIAG2)
	END IF
 24	FORMAT(/1X,'CELL DISTANCE CHECK (DIAGONAL LENGTH) IS ',F5.2)
C
C	NOW READ IN COORDINATES AND THEN OPERATE ON THEM.
C	LAST RECORD IS 'END'
	I=0
 100	I=I+1
	READ (IUNIT,105,END=110) ATNAME(I),(X(J,I),J=1,3)
	GO TO 100
 105	FORMAT(A,5X,3F10.5)
 110	CONTINUE
	IATOM=I-2
C
C	  HOW MANY NEIGHBORING CELLS TO GENERATE?
C	  IF WE'RE DOING A DISTANCE CHECK, WE NEED TO USE 5X5 ARRAY
C	  IN CASE FRAC. COORDS WERE GIVEN OUTSIDE OF -1/+1 RANGE,
C	  AS THEY SOMETIMES ARE.
C
	IF(ALL27)THEN
		IF(DCHECK)THEN
			NNBRS=7
		ELSE
			NNBRS=3
		END IF
	ELSE
		NNBRS=1
	END IF
	N=0
	NGEN=0
	MWRITN=0
	NWRIT=0
	WRITEM=.NOT.DCHECK	!DON'T OUTPUT ASYMM. UNIT FOR OPT=3
C
	DO IX=1,NNBRS
	UTR(1)=NBC(IX)
	DO IY=1,NNBRS
	UTR(2)=NBC(IY)
	DO IZ=1,NNBRS
	UTR(3)=NBC(IZ)
C
	N=N+1
C
C	  COMPUTE CENTROID OF DISTANCE TEST MOLECULE,
C	  BUT ONLY AFTER FIRST PASS (AFTER N=1=UNIT CELL).
C
	IF(DCHECK .AND. N.EQ.2)THEN
		XMIN=UX(1,1)
		YMIN=UX(2,1)
		ZMIN=UX(3,1)
		XMAX=XMIN
		YMAX=YMIN
		ZMAX=ZMIN
		DO K=2,IATOM
			XMIN=MIN(XMIN,UX(1,K))
			YMIN=MIN(YMIN,UX(2,K))
			ZMIN=MIN(ZMIN,UX(3,K))
			XMIN=MIN(XMIN,UX(1,K))
			YMIN=MIN(YMIN,UX(2,K))
			ZMIN=MIN(ZMIN,UX(3,K))
		END DO
		XCT=(XMAX+XMIN)/2.
		YCT=(YMAX+YMIN)/2.
		ZCT=(ZMAX+ZMIN)/2.
	END IF
C
C	  COMPUTE DISTANCE OF EDEGS OF THIS CELL FROM 000.
C	  IF TOO FAR, SKIP THIS CELL COMPLETELY.
C
	IF( DCHECK .AND. (N.GT.1) )THEN
		CALL MATMPY(TMAT,UTR,CCORN,TTR)
		DIST2=(CCORN(1)-XCT)**2 + (CCORN(2)-YCT)**2 + (CCORN(3)-ZCT)**2
		IF(DIST2.GT.DIAG2)THEN
			SKIP=.TRUE.
		ELSE
			SKIP=.FALSE.
		END IF
	ELSE
		SKIP=.FALSE.
	END IF
C
	CWRIT=.FALSE.
	IF(.NOT.SKIP)THEN
	NGEN=NGEN+1
	DO J=1,NSYMM
		DO I=1,IATOM
			CALL MATMPY(SYMM(1,1,J),X(1,I),AX,STR(1,J))
C
			IF(BOUNDS)THEN
				DO K=1,3
				IF( AX(K).LT. 0. )AX(K)=AX(K)+1.
				IF( AX(K).GT. 1. )AX(K)=AX(K)-1.
				END DO
			END IF
C
			CALL MATMPY(UNIT,AX,BX,UTR)
			CALL MATMPY(TMAT,BX,CX(1,I),TTR)
C
C		SAVE ASYMM. UNIT CART. COORDS. FOR LATER DISTANCE CHECK.
C
			IF( (N.EQ.1) .AND. (J.EQ.1) )THEN
				UX(1,I)=CX(1,I)
				UX(2,I)=CX(2,I)
				UX(3,I)=CX(3,I)
			ELSE
				IF(DCHECK)THEN
					SKIPATOM = HEAVY .AND. (ATNAME(I)(1:1).EQ.'H')
					JJ=0
					DO WHILE ((WRITEM.NE. .TRUE.).AND.
     &						  (JJ.LE.IATOM) .AND. .NOT.SKIPATOM )
					JJ=JJ+1
					NN=(JJ-1)/IATOM
					II=JJ-NN*IATOM
					SKIPATOM = HEAVY .AND. (ATNAME(II)(1:1).EQ.'H')
					IF(.NOT.SKIPATOM)THEN
						D2=( (UX(1,JJ) - CX(1,I))**2 +
     &						     (UX(2,JJ) - CX(2,I))**2 +
     &						     (UX(3,JJ) - CX(3,I))**2 )
						IF(D2.LT.D2MAX)THEN
							WRITEM=.TRUE.
							WRITE(6,202)(IFIX(UTR(K)),K=1,3),
     &							J,I,SQRT(D2),NN+1,II,SQRT(DIST2)
						END IF
					END IF
					END DO
				ELSE
					WRITEM=.TRUE.
				END IF
			END IF   !DON'T EVEN BOTHER DISTANCE CHECK FOR ASYMM. UNIT
 202	FORMAT(/1X,'CELL ',3I2,'  SYMMETRY ELEMENT ',I2,'   ATOM #',I3,
     &		' IS ',F5.2,' ANGSTROMS FROM'/
     &	        1X, 'CELL  0 0 0  SYMMETRY ELEMENT ',I2,'   ATOM #',I3,
     &		/1X,'CELL DISTANCE = ',F5.2)
C
		END DO  ! LOOP OF #ATOMS IN ASYMMETRIC UNIT
C
C	   NOW WRITE OUT THIS ASYMMETRIC UNIT IF NEEDED.
C	   CWRIT REMEMBERS IF THIS CELL # HAS BEEN WRITTEN OUT YET.
C
		IF(WRITEM)THEN
			IF(.NOT.CWRIT)THEN
				CWRIT=.TRUE.
				NWRIT=NWRIT+1
			END IF
C			WRITE (6,4) AX,BX,CX
			MWRITN=MWRITN+1
			WRITE (OUNIT,44) (ATNAME(I)(1:4),ATOMID(NWRIT,J),(CX(K,I),K=1,3),
     &					I=1,IATOM)
			CALL MOTIFPIX (CX,IATOM)
			WRITEM=.FALSE.
		END IF
	END DO   ! LOOP OVER # OF SYMMETRY ELEMENTS
C
	END IF  !IF DISTANCE CHECK MADE US SKIP OVER THIS CELL
	END DO  !X TRANSLATION OF UNIT CELLS
	END DO  !Y     "
	END DO  !Z     "
 4	FORMAT (//1X,3F10.5,/1X,3F10.5,/1X,3F10.5)
 44	FORMAT (A,A,3X,3F10.5)
C
	WRITE(OUNIT,54)
 54	FORMAT('END')
 999	WRITE(6,55)N,NGEN,NWRIT,MWRITN
 55	FORMAT(/1X,I3,' CELLS WERE CONSIDERED.'
     &	       /1X,I3,' CELLS HAD COORDINATES GENERATED.'
     &	       /1X,I3,' CELLS HAD (SOME) COORDINATES OUTPUT.'
     &	       /1X,I3,' TOTAL MOLECULES WRITTEN.')
	STOP
C
	END

C==============================================================================
	SUBROUTINE MATMPY(ROT,AX,BX,TR)
C
	DIMENSION ROT(3,3),AX(3),BX(3),TR(3)
C
	DO I=1,3
		BX(I)=TR(I)
		DO J=1,3
			BX(I)=BX(I)+ROT(I,J)*AX(J)
		END DO
	END DO
C
	RETURN
C
	END
C

C==============================================================================
	CHARACTER*2 FUNCTION ATOMID(N,J)
	CHARACTER*1 NCELL(63)
	PARAMETER MAXSYM=20
C
C	  MAKES AN ATOM IDENTIFIER FROM N (UNIT CELL NUMBER) AND
C					J (SYMMETRY OPERATION NUMBER)
C
	NCELL(1)=' '
	DO I=1,26
	NCELL(I+1)=CHAR(I-1+ICHAR('A'))
	END DO
	DO I=1,26
	NCELL(I+27)=CHAR(I-1+ICHAR('a'))
	END DO
	DO I=1,10
	NCELL(I+53)=CHAR(I-1+ICHAR('0'))
	END DO
C
C	I=(N-1)*MAXSYM + J
C	K1=(I-1)/63+1
C	K2=I-63*(K1-1)
C
C	ATOMID=NCELL(K1)//NCELL(K2)
C
	IF(N.GT.63)THEN
		WRITE(6,*)'TOO MANY CELLS GENERATED!!!'
		ATOMID='??'
	ELSE
		ATOMID=NCELL(N)//CHAR( J + ICHAR('0') )
	END IF
C
	RETURN
	END
C

C==============================================================================
	SUBROUTINE TJMATRIX(CELL,VOL,T)
C. . . . SETS UP THE ORTHOGONALIZATION MATRIX TO GO FROM
C. . . . FRACTIONAL COORDINATES IN x,y,z (CRYSTAL SPACE)
C. . . . TO AN ORTHOGONAL FRAME OF REFERENCE X,Y,Z 
C. . . . NEW AXIS WILL BE:
C        X' ALONG A-AXIS, Y' IN THE (A,B) PLANE AND Z' ALONG C*.
C
        DIMENSION CELL(6),T(3,3)
C
        A = CELL(1)
        B = CELL(2)
        C = CELL(3)
        COSA= CELL(4)
        COSB= CELL(5)
        COSG= CELL(6)
        SINA = SQRT(1.-COSA*COSA)
        SINB = SQRT(1.-COSB*COSB)
        SING = SQRT(1.-COSG*COSG)
        T(1,1) = A
        T(1,2) = B*COSG
        T(1,3) = C*COSB
        T(2,1) = 0.
        T(2,2) = B*SING
        T(2,3)=C*(COSA-COSB*COSG)/SING
        T(3,1) = 0.
        T(3,2) = 0.
        T(3,3) = VOL/(A*B*SING)
        RETURN
        END
C

C==============================================================================
	REAL FUNCTION VOL(CELL)
C
C . . . GIVEN CELL COMPUTES VOLUME.
C
        DIMENSION CELL(6)
C . . . .
        COSA= CELL(4)
        COSB= CELL(5)
        COSG= CELL(6)
        SINA= SQRT(1.-COSA*COSA)
        SINB= SQRT(1.-COSB*COSB)
        SING= SQRT(1.-COSG*COSG)
        A1 = CELL(1)
        A2 = CELL(2)
        A3 = CELL(3)
        TEMP =(1.-COSA*COSA-COSB*COSB-COSG*COSG+2.*COSA*COSB*COSG)
        VOL =A1*A2*A3*SQRT(TEMP)
C
	RETURN
	END
C

C==============================================================================
	FUNCTION TRIM(STRING)
C
	IMPLICIT INTEGER (A-Z)
	CHARACTER STRING*(*)
C
	DO 10 I=LEN(STRING),1,-1
	IF(STRING(I:I).NE.' ')THEN
		TRIM=I
		RETURN
	END IF
 10	CONTINUE
C
	TRIM=0
	RETURN
	END	

C==============================================================================
	SUBROUTINE MATVEC1(C,A,B)
C . . ROUTINE TO PRE-MULTIPLY A VECTOR B(3), BY A MATRIX (3,4)
C . . TO YIELD A VECTOR C(3)
C . .      C(I) = A(I,K)*B(K)
C . .
		DIMENSION A(3,4),B(3),C(3)
C . .
               
		DO I=1,3
		  C(I) = A(I,4)
		  DO J=1,3
			C(I) = C(I) + A(I,J)* B(J)
		  END DO
		END DO
		RETURN
		END

C==============================================================================
		SUBROUTINE MATVEC2(C,A,B)
C . . ROUTINE TO PRE-MULTIPLY A VECTOR B(3), BY A MATRIX (3,3)
C . . TO YIELD A VECTOR C(3); NO TRANSLATIONS
C . .      C(I) = A(I,K)*B(K)
C . .
		DIMENSION A(3,3),B(3),C(3)
C . .
               
		DO I=1,3
		  C(I) = 0.0
		  DO J=1,3
			C(I) = C(I) + A(I,J)* B(J)
		  END DO
		END DO
		RETURN
		END

C==============================================================================
                SUBROUTINE MATPRO(CM,AM,BM)
C . . ROUTINE TO MUTIPLY TWO-SQUARE MATRICES 3 BY 3 IN THE
C . . ORDER CM(I,J) = AM(I,K)*BM(K,J)
C . .
                DIMENSION CM(3,3),AM(3,3),BM(3,3)
C . .
		DO I=1,3
		   DO J=1,3
		     CM(I,J) = 0.0
			DO K=1,3
		          CM(I,J) = CM(I,J) + AM(I,K)*BM(K,J)
			END DO
                   END DO
                END DO
		RETURN
		END

C==============================================================================
      SUBROUTINE READPDB(ONLYCA,OKAY,IPDB,JCRY,ROT,TM,TMI,XOR,NOR)
C . . INITIAL CODE FOR THIS ROUTINE WAS OBTAINED FROM THE
C . . PROTEIN DISPLAY SYSTEM (PDS) OF ABBOTT LABORATORIES
C . . DECEMBER 15. 1986. LATER MODIFIED BY C. ABAD-ZAPATERO
C . . TO FIT THE NEEDS OF THE  TABLES  PROGRAM.

C	SUBROUTINE PF_BKNIN (OKAY, MOL, MOLNAME, UNIT, NONEW)
C	*****************************************************
C       ORIGINAL ROUTINE NAME IN PDS PF_BKNIN
C
C	READ IN A BROOKHAVEN FILE
C
C	OKAY		- TRUE IF SUCCESSFUL
C	IPDB		- UNIT TO READ (ALREADY OPENED)
C	NONEW		- NO NEW LIBRARY RESIDUES
C
	IMPLICIT NONE
	INCLUDE 'PARAMS.INC'

C . .
	LOGICAL		OKAY,NONEW,ONLYCA
	INTEGER		MOL,NUIN,NUAU,NID,NCA,NMOD,IPDB,NOR
        INTEGER         M,IOR,K,JCRY 
C
C . . THE MATRICES ROT,TM,TMI AND XOR ARE TRANSFERRED FROM THE MAIN PROGRAM:
       REAL*4 ROT(3,4),TM(3,3),TMI(3,3),XOR(3,4)
C . . ARRAYS VI,VO,RF,RO ARE WORKING ARRAYS FOR INPUT AND OUTPUT OF COORDINATES:
       REAL*4 VI(3),VO(3),RF(3),RO(3)
C . .
       CHARACTER*4 WHICH
C . .
	REAL		AXYZS(3,MAX_ATMRES)
	INTEGER		MAXCONECT
	PARAMETER	(MAXCONECT = MAX_CROSSES * 4)
	INTEGER		ISTAT, IA, NHEADER, NCOMPND, NSOURCE, NAUTHOR
	INTEGER		NSSBOND, NATOM, NTER, NCONECT, NCROSS, NCHN
	INTEGER		NAOKAY, NADONE, ICOMP, IRES, ICHN, IABASE, ITMP
	INTEGER		J, RESNUM, NATMRES
	INTEGER*2	ASERIALS(MAX_ATMRES)
	INTEGER		ATMS(MAX_ATMRES), ALLATMS(MAX_ATOMS)
	INTEGER		SSCIDX(MAXCONECT)
	INTEGER		NCTOT
	INTEGER		CCHAIN(2,MAXCONECT)
	INTEGER		CRESNUM(2,MAXCONECT)
	INTEGER		CONNEI(4,MAXCONECT)
	INTEGER		CONATM(MAXCONECT)
	LOGICAL		OKAY2
	LOGICAL		NMD, OD
	LOGICAL		NEWCHN
	CHARACTER*(LEN_LINE)	ERRLINE, INLINE
	CHARACTER*128	HDR, CMPD, SRCE, AUTH
	CHARACTER*16	TSTR, LTYPE
	CHARACTER*16	MODNAME
	CHARACTER*(LEN_LINE)	FNAME
	CHARACTER*16	FM, SEQ
	CHARACTER*(LEN_FIELD)	CHNNAM(MAX_CHAINS)
	CHARACTER*1	CHNID
	CHARACTER*4	RESNAM
	CHARACTER*1	RESLET
	CHARACTER*4	ANAMES(MAX_ATMRES)
	CHARACTER*1	SSCHAIN(2,MAXCONECT)
	CHARACTER*1	CRESLET(2,MAXCONECT)
	CHARACTER*4	CATMNAM(2,MAXCONECT)
C
C	*** BROOKHAVEN ATOM,HETATOM RECORD (INCOMPLETE)
	INTEGER		BSERIAL
	CHARACTER*4	BNAME,FIN
	CHARACTER*1	BALTLOC
	CHARACTER*4	BRESNAM
	CHARACTER*1	BCHNID,IATN
	INTEGER		BRESNUM
	CHARACTER*1	BRESLET
	REAL		BXYZ(3)
        DATA IATN /'C'/
        DATA FIN /'END '/
C
C	*** INITIALIZE AND ERROR CHECKS
	OKAY = .FALSE.
	INQUIRE (IPDB, ERR=900, IOSTAT=ISTAT,
	1	FORM=FM, NAME=FNAME, NAMED=NMD, OPENED=OD, SEQUENTIAL=SEQ)
	IF (.NOT. OD) THEN
		ERRLINE = 'File IPDB not open: unit 3 '
		GOTO 999
	END IF
	IF (.NOT. NMD) THEN
		ERRLINE = 'File IPDB opened without name: unit 3 '
		GOTO 999
	END IF
	IF (SEQ .NE. 'YES') THEN
		ERRLINE = 'IPDB not opened for sequential access '
		GOTO 999
	END IF
	IF (FM .NE. 'FORMATTED') THEN
		ERRLINE = 'IPDB not opended for FORMATTED access '
		GOTO 999
	END IF
C
C	*** SETUP FOR READING
	NHEADER = 0
	NCOMPND = 0
	NSOURCE = 0
	NAUTHOR = 0
	NSSBOND = 0
	NATOM = 0
	NTER = 0
	NCONECT = 0
	NCROSS = 0
C
C . . . .ZERO OUT A FEW VARIABLES
        NUIN = 0
        NUAU = 0
        NID = 0
        NCA = 0
C . .EVERY NMOD(THG) ATOM WILL BE PRINTED
        NMOD=10  
C . .
	NCHN = 0
	NEWCHN = .TRUE.
	CHNID = CHR_BLANK
	RESNUM = WLD_RESN
	RESLET = CHR_WILD
	NATMRES = 0
	NCTOT = 0
	NAOKAY = 0
	NADONE = 0
C
C	*** READ IN THE FILE
100	CONTINUE
	    READ (IPDB, 101, END=199, ERR=910, IOSTAT=ISTAT) INLINE
101	    FORMAT (A)
C
C	    *** PARSE THE LINE
	    LTYPE = INLINE(1:6)
	    IF (LTYPE .EQ. 'HEADER') THEN
		NHEADER = NHEADER + 1
		IF (NHEADER .EQ. 1) HDR = INLINE(11:60)
	    ELSE IF (LTYPE .EQ. 'COMPND') THEN
		NCOMPND = NCOMPND + 1
		IF (NCOMPND .EQ. 1) CMPD = INLINE(11:70)
	    ELSE IF (LTYPE .EQ. 'SOURCE') THEN
		NSOURCE = NSOURCE + 1
		IF (NSOURCE .EQ. 1) SRCE = INLINE(11:70)
	    ELSE IF (LTYPE .EQ. 'AUTHOR') THEN
		NAUTHOR = NAUTHOR + 1
		IF (NAUTHOR .EQ. 1) AUTH = INLINE(11:70)
	    ELSE IF (LTYPE .EQ. 'SSBOND') THEN
		NSSBOND = NSSBOND + 1
		NCTOT = NCTOT + 1
		IF (NCTOT .GT. MAXCONECT) THEN
		    ERRLINE = 'Too many crosslinks: '//INLINE
		    GOTO 999
		END IF
		READ (INLINE, 110, ERR=920, IOSTAT=ISTAT)
	1	    SSCHAIN(1,NSSBOND),CRESNUM(1,NCTOT),CRESLET(1,NCTOT),
	2	    SSCHAIN(2,NSSBOND),CRESNUM(2,NCTOT),CRESLET(2,NCTOT)
110		FORMAT (6X,1X,3X,1X,3X,1X,A1,1X,I4,A1,3X,3X,1X,A1,1X,I4,A1)
		CCHAIN(1,NCTOT) = -1
		CCHAIN(2,NCTOT) = -1
		CATMNAM(1,NCTOT) = 'SG'
		CATMNAM(2,NCTOT) = 'SG'
		SSCIDX(NSSBOND) = NCTOT
C . . .
C . . . THIS IS THE IMPORTANT SWITCH FOR PACKING STUDIES; READ ATOMS
C . . . SELECT ACCORDING TO THE ONLYCA LOGICAL VARIBLE AND THEN
C . . . WRITE OUT ON SCRATCH FILE JCRY FOR LATER USE. .
	    ELSE IF ((LTYPE .EQ. 'ATOM') .OR. (LTYPE .EQ. 'HETATM')) THEN
		NUIN = NUIN + 1
		READ (INLINE, 111, ERR=920, IOSTAT=ISTAT)
	1		BSERIAL, BNAME, BALTLOC, BRESNAM,
	1		BCHNID, BRESNUM, BRESLET, BXYZ
111		FORMAT (6X,I5,1X,A4,A1,A4,A1,I4,A1,3X,3F8.3)
C . . . BASIC INFORMATION ABOUT THE ATOM READ; CHECK FOR CA*S AND
C . . . WRITE OUT:
C . . 
        IF (BNAME.NE.' CA ' ) GO TO 198
C .. APPLY ROTATION AND TRANSLATION IN ORTHOGONAL SET OF AXIS
C .. VECTOR POST-MULTIPLIES MATRIX   XX= ROT(I,J)*X
        NCA = NCA + 1
        VI(1) = BXYZ(1)
        VI(2) = BXYZ(2)
        VI(3) = BXYZ(3)
        CALL MATVEC1(VO,ROT,VI)
C . . . 
C . . . CHANGE COORDINATES TO FRACTIONAL COORDINATES IN
C ..... OBLIQUE SYSTEM . . . . . . . . .
C . . . . MATRIX TMI (ABOVE) PRE-MULTIPLIES VECTOR 
C . . . . ALGEBRA FROM  NOTES ON MOLECULAR
C . . . . REPLACEMENT FROM J.E. JOHNSON AT PURDUE UNIVERS.
C . . . .  x(fractional) = TMI*X (angstroms orthogonal)
        DO M=1,3
          VI(M) = VO(M)
          RF(M) =0.
          RO(M) = 0.
        END DO
          CALL MATVEC2(RF,TMI,VI)
C . . .
C . . .         
C . . . RF(3) CONTAINS NOW THE FRACTIONAL COORDINATES
C . . . LOOP TO TAKE INTO ACCOUNT THE DIFFERENT BRAVAIS LATTICES
        DO IOR=1,NOR
          NUAU= NUAU + 1
          NID = NUAU +99 
             DO K=1,3
               RF(K) = XOR(K,IOR) + RF(K)
             END DO
C . . WRITE OUT THE TEMPORARY FILE SPG1.CRY CONTAINING
C . . FRACTIONAL COORDINATES OF ATOMS  
          WRITE(JCRY,800)IATN,NID,RF
          IF(MOD(NUIN,NMOD).EQ.0)
     .    WRITE(6,25) BXYZ,VO,RF,BNAME
 25       FORMAT(2X,9F10.5,1X,A4)
 800      FORMAT(A1,I3,6X,3F10.5)
        END DO 

	    ELSE IF (LTYPE .EQ. 'USER') THEN
		LTYPE = INLINE(8:13)
		IF (LTYPE .EQ. 'CROSS') THEN
		    NCROSS = NCROSS + 1
		    NCTOT = NCTOT + 1
		    IF (NCTOT .GT. MAX_CROSSES) THEN
			ERRLINE = 'Too many crosslinks: '//INLINE
			GOTO 999
		    END IF
		    READ (INLINE, 113, ERR=920, IOSTAT=ISTAT)
	1		CCHAIN(1,NCTOT), CRESNUM(1,NCTOT),
	1		CRESLET(1,NCTOT), CATMNAM(1,NCTOT),
	1		CCHAIN(2,NCTOT), CRESNUM(2,NCTOT),
	1		CRESLET(2,NCTOT), CATMNAM(2,NCTOT)
113		    FORMAT (6X,1X,6X,1X,I2,1X,I4,A1,1X,A4,3X,I2,1X,I4,A1,1X,A4)
		ELSE
C		    *** UNRECOGNIZED USER LINE TYPE
		END IF
	    ELSE
C		*** UNRECOGNIZED BROOKHAVEN LINE TYPE
	    END IF
C
198	    CONTINUE
	    GOTO 100
199	CONTINUE
C
C	*** ALL DONE
	OKAY = .TRUE.
        IF (OKAY) THEN
         WRITE(6,30) NUIN,NCA
 30      FORMAT(5X, 'NUMBER OF ATOMS READ IN',2X,I10,/
     .   ,5X, 'NUMBER OF CA ATOMS WRITTEN OUT',2X,I10)
        END IF        
        WRITE(JCRY,6990) FIN
 6990   FORMAT(A)
        CLOSE (UNIT=JCRY)
	RETURN
C
C
C	*** ERROR TRAPS
900	CONTINUE
	ERRLINE = 'Attempted to inquire about unit 3: IPDB file'
	GOTO 999
C
910	CONTINUE
	ERRLINE = 'Reading line in unit 3: IPDB file'
	GOTO 999
C
920	CONTINUE
	ERRLINE = 'Decoding line in unit 3: IPDB  file'
	GOTO 999
C
999	CONTINUE
        WRITE(6,3000)ERRLINE
3000    FORMAT(A) 
	RETURN
C
	END
C . .
C . .
C==============================================================================
        SUBROUTINE OUTCELLPIX
     .  (IO,DAT,TIM,PROG,VERS,IAB,NA,NB)
C . .  WRITES THE CELL.PIX FILE ON UNIT IO
C . .  AND RETURNS
C . .  BASED ON TRANSFERRED INFORMATION IN CALLING ROUTINE
C . . 
C . .  THE IMPORTANT VALUES: NUMBER OF VERTICES (ATOMS): NA 
C . .                      : NUMBER OF LINES (BONDS): NB
C . .
      CHARACTER*1 C
      CHARACTER*4 BONDS
      INTEGER*2 IBONDS
C . .
C . . SET COMMON BLOCKS FOR THE CELL FILE:
C . . 
C . . 
      COMMON /CELL/V000(3,8),BONDS(2,12),IBONDS(2,12),C,WG1,BONTY1     
C . .
        CHARACTER*9 PROG,VERS
        CHARACTER*8 DAT,TIM
        INTEGER*2 IX,IY,IZ,IW
C
C . . THE VECTORS TO MAKE THE PARALLELEPIPED.
C . . THE CONNECTIVITY IS STORED IN BONDS ARRAY.
C . . THE FACTOR OF 1000. IS TO SCALE ANGSTROMS TO GRAMPS' 16BIT
C . . INTEGER SCALE, WHERE MAXIMUM COORDINATE IS 32767.  IF YOUR
C . . COORDINATES ARE GREATER THAN 32.767 A, THEN THIS WILL NOT
C . . WORK, SO CHANGE THE 1000. FACTOR TO WHATEVER YOU'D LIKE.
         WRITE(IO,*) NB*2,0,4
         DO IB=1,NB
         DO K=1,2
         IX=V000(1,IBONDS(K,IB)+1)*1000.
         IY=V000(2,IBONDS(K,IB)+1)*1000.
         IZ=V000(3,IBONDS(K,IB)+1)*1000.
         WRITE(IO,*) IX,IY,IZ
         END DO
         END DO
C . . ALL DONE! 
        RETURN
        END
C . .

C==============================================================================
        SUBROUTINE OUTOPSPIX
     .  (IO,DAT,TIM,PROG,VERS,IAB,NSYP,NSYV)
C . .  WRITES THE THE SYMOPS.PIX TYPE FILE ON UNIT IO:
C . .  BASED ON TRANSFERRED INFORMATION IN CALLING ROUTINE
C . .  AND COMMON BLOCK SYMOPS.
C . .  THE IMPORTANT VALUES: NUMBER OF SYMMETRY POINTS (ATOMS): NSYP 
C . .                      : NUMBER OF SYMMETRY VECTORS NSYV(BONDS):NSYV
C . .
      PARAMETER MAXATOM=2000
      PARAMETER MAXSYM =32
      PARAMETER MAXOPS =200
C . .
        CHARACTER*1 IATID
        INTEGER*2 IX,IY,IZ,IW
C . .
C . . SET COMMON BLOCKS FOR THE SYMOPS FILE : MAXOPS = 200
C . . 
      COMMON /SYMOPS/SVEC(3,MAXOPS),IATID(MAXOPS),WG2,BONTY2
C . .
        CHARACTER*9 PROG,VERS
        CHARACTER*8 DAT,TIM
C . .
C . . . WRITE THE VECTORS DEFINING THE DIRECTIONS OF THE
C . . . SYMMETRY ELEMENTS IN THE CELL
          WRITE(IO,*)NSYV*2,0,4
          DO NV=1,NSYV
           I1 = 1 + (NV-1)*2
           I2 = I1 + 1
            IX=SVEC(1,I1)*1000.
            IY=SVEC(2,I1)*1000.
            IZ=SVEC(3,I1)*1000.
            WRITE(IO,*)IX,IY,IZ
            IX=SVEC(1,I2)*1000.
            IY=SVEC(2,I2)*1000.
            IZ=SVEC(3,I2)*1000.
            WRITE(IO,*)IX,IY,IZ
          END DO
C . . ALL DONE! 
        RETURN
        END
C . .
	SUBROUTINE MOTIFPIX (CX,NATOMS)
	REAL CX(3,NATOMS)
	INTEGER NATOMS
C
C	  ROUTINE TO WRITE OUT A PIX FILE OF A MOLECULE OF NATOMS
C	  WHOSE COORDINATES ARE IN CX.
C
	PARAMETER PUNIT=13
	PARAMETER MAXBONDS=100
	INTEGER IJ(MAXBONDS*2)
	REAL DIST2
	REAL BONDED /1.7/  !ANGSTROMS
C .
C .
C . OPEN FILE SPG2.PIX ==> CRYSTAL PACKING IN PIX FORMAT
C .
        OPEN(UNIT=PUNIT,FILE='SPG2.PIX',STATUS='NEW',SHARED)

	IBOND=1
	DO I=1,NATOMS-1
	DO J=I+1,NATOMS
	IF( ABS(CX(1,I)-CX(1,J)) .LE. BONDED .OR.
     &	    ABS(CX(2,I)-CX(2,J)) .LE. BONDED .OR.
     &	    ABS(CX(3,I)-CX(3,J)) .LE. BONDED )THEN
		DIST2=0.
		DO K=1,3
		DIST2=DIST2+(CX(K,I)-CX(K,J))**2
		END DO
		IF(DIST2.LE.BONDED**2)THEN
			IJ(IBOND)=I
			IJ(IBOND+1)=J
			IBOND=IBOND+2
		END IF
	END IF
	END DO
	END DO
C
	WRITE(PUNIT,*)IBOND-1,0,4
	DO K=1,IBOND-1,2
	WRITE(PUNIT,*)	( INT(CX(L,IJ(K  ))*1000.),L=1,3 ),IJ(K)
	WRITE(PUNIT,*)	( INT(CX(L,IJ(K+1))*1000.),L=1,3 ),IJ(K+1)
	END DO
C
	RETURN
	END
