      PROGRAM TRACOR
**************************   File: TRACOR.NEW: 21 July '97    TRACOR.FOR
*                                            Last update:   11 Nov. 1998
*                                            Source: ROG 1986 , SGG 1987
* After 'CSUBFILE FOUR2' follows the second part of TRACOR = Fourier pgm
************************************************************************
C wensen: -  FOUR groter in Tracor    - FOM voor BOTS !
 
*TRACOR LOG of recent modifications (last on top
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      EQUIVALENCE (ICRYS, IFILE(3))
      EQUIVALENCE (IBINFF, IFILE(16))
      EQUIVALENCE (LIS1, IFILE(7)), (LIS2, IFILE(8))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      DIMENSION SEDEL2(50,2)
      CALL KEPROG('TRACOR')
      WRITE (LIS2, FMT = '('' Last TRACOR update: 11 Nov. 1999'')')
      CALL RDCRYS (ICRYS)
      IF (IPOLA .EQ. 7) CALL KERROR
     *      ('TRACOR not applicable to space group P1' , 0, 'TRACIN')
      CALL FCALIN
      CALL TRACIN
      CALL COTRA1
      IF (IPOLA .LE. 0) THEN
         CALL D3FOUR
         CALL FILCLO (IBINFF, 'DELETE')
      ELSE
         CALL D12FOU
         ENDIF
      CALL COTRA3
      CALL R2TRAC (SEDEL2)
      CALL SHIFTR (SEDEL2)
      CALL KEPROX
      WRITE (LIS2, FMT='(/'' Test: DDOKA exit MAIN SUBPROGRAM ''/)')
      STOP 99
      END
      SUBROUTINE FCALIN
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      EQUIVALENCE (IDDL, IFILE(1)), (IATOMS, IFILE(2))
      EQUIVALENCE (IDDS, IFILE(1)), (ICRYS,IFILE(3))
      EQUIVALENCE (ICON, IFILE(4)), (LIS1, IFILE(7)), (LIS2, IFILE(8))
      EQUIVALENCE (IBINFO, IFILE(11)), (IBINFC, IFILE(12))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     *               WAVE,     CELALL(10),  AMOLW,      ZET,
     *               NELEC,    F000,        ABSMU,      ICENT,
     *               ILATT,    ISYST,       ILAUE,      IMULT,
     *               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     *         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     *         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      COMMON /CRYSB/ SPGR,     WAVEAT,      CELATY(10)
      CHARACTER      SPGR *16, WAVEAT *2,   CELATY *2
      COMMON /FCALCA/ BP,       BR,       SCALE,    HKLMAX(3), STLMAX,
     *                IZTYPE(10), CELPAR(10), PSQ,  P1SQ,     ITRS(24),
     *        AMULT,  ASYMM,    ALATT,    ASYMCL,   NSYMC,    ASYMC,
     *                HKLX(3,24), IDHKL(24), HCODE, FOBS,     SIG,
     *                STL,      STL2,     ISS,      ENORM,
     *                FP,       PHIP,     FAP,      FBP,      EPSIL,
     *                EPSIL2,   SF2,      SF2P,     FPEXP(2,24)
      DIMENSION FITFO(3), FITFC2(51)
      EQUIVALENCE (HCODE, FITFO(1)), (EPSIL2, FITFC2(1))
      COMMON /COTRA/ MULS(48,48), INDX(3,3), IDIMF, TT(3,48),
     *               EFST(48), FHST(48), PHST(48), CELLN(6), MAXHKL(3),
     *  P111, P77, EMIN, EO2AV, SCX, BOVX, DAMPX, SMM, IOMAPS, SMAX,
     *  IRR(3,3,48), IORG, BBB, D2R, R2D,
     *  ATEM1, ATEM2, ATEM3, ATEM4, NTEM5, ATEM6
      PARAMETER (MAXAT = 993)
      PARAMETER (MAXBUF = 198)
      COMMON /  / SICO(12500), FF(500,10),  EXPBP(500),   EXPBR(500),
     *            SUMF2(500),  SUMF2P(500), SFAC(13,10),
     *            ATXYZ(10,MAXAT), IZAT(MAXAT), ITAT(MAXAT), NAT,
     *            BUFFO(MAXBUF),  BUFFC(MAXBUF), BUFBUF(MAXBUF)
      DIMENSION BLACOM(42000)
      EQUIVALENCE (BLACOM(1), SICO(1))
      COMMON /ATNAMA/ ATNAME(MAXAT)
      CHARACTER * 6   ATNAME
      PARAMETER (LCMAX = 15)
      CHARACTER * 6 LCONDA(LCMAX)
      DATA LCONDA / 'TRACOR', 'EMIN'  , 'SCSG'  , 'BHSG'  , 'DAMP'  ,
     *              'SCALE' , 'BBB'   , 'STLMAX', 'SMM'   , 'MAXHKL',
     *              'PRIMAP', 'PSQMAX', 'XXXXXX', 'WILSON', 'PRINT' /
      WRITE (LIS2, 101)
  101 FORMAT (' Structure factor calculation for TRACOR:'
     *      / ' EXPAND data to P1 symmetry (or centered equivalent)')
      CALL RDCRYS (ICRYS)
      DO 110 I=1,NTYPE
  110 BUFFO(I) = CELALL(I) / ZET
      I = NTYPE
      J = NINT(ZET)
      WRITE (LIS2, 114) J, (CELATY(K), BUFFO(K), K=1,I)
  114 FORMAT (' Z:', I3 / ' FORMUL:', 6(2X,A2,F6.1) /
     *                           ( 8X, 6(2X,A2,F6.1)))
      CALL LOGRD (IDDL, 'MERBSC', KLOG)
      IF (KLOG .LT. 0 .OR. LIT(2) .NE. 'SCALE') THEN
         CHOUT=' MERBIN SCALE not found on DDLOG file: what happened?'
         CALL SHOUT
         CALL KERROR ('You cleared too much, probably!', 114, 'FCALIN')
         ENDIF
      SCALE = FNUM(2)
      BOV = FNUM(3)
      BP = BOV
      BR = BOV
      WRITE (LIS2, 116) SCALE, BOV
  116 FORMAT (' Data from DDLOG: Scale =', F9.5, ' Bov:', F6.3)
      CALL FILCLO (IDDL, 'KEEP')
      CALL BINIFF (1, IBINFO, 'BINFO', FITFO, NITFO, BUFFO, KENDFO)
      STLMAX = BUFFO(6)
      CALL KERNAB (BUFFO(7), HKLMAX, 3)
      CALL KERF2I (HKLMAX, MAXHKL, 3)
      EMIN = 0.0
      SCX = SCALE
      BOVX = BOV
      DAMPX = BP
      SMM  = 0.00001
      IOMAPS = 0
      P77 = 1.
  120 CALL RDCOND (ICON, LCONDA, LCMAX, KEND)
      GOTO (120, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 120, 120,120), KEND
      IF (KEND.LE.0) GOTO 140
  2   EMIN = FNUM(1)
      WRITE (LIS1, FMT = '('' Input value for Emin: '', F10.2)') EMIN
      GOTO 120
  3   SCSG = FNUM(1)
      IF(SCSG .LT. 0.001) SCSG = 1.0000
      SCX = SCALE * SCSG
      WRITE (LIS1, FMT = '('' SCALE '', F10.5,'' = SCALE '', F10.5,
     * '' multiplied by extra scale: '', F10.5)') SCX, SCALE, SCSG
      GOTO 120
  4   BHSG=FNUM(1)
      BOVX = BOV + BHSG
      WRITE (LIS1, FMT = '('' Input: additional Bp'', F10.5 )') BHSG
      WRITE (LIS1, FMT = '('' Now to be used:  Bov'', F10.5 )') BOVX
      GOTO 120
  5   DAMP = 0.5 * FNUM(1)
      DAMPX = DAMPX - DAMP
      WRITE (LIS1, FMT = '('' Supplied: damping parameter '', F10.5)')
     * FNUM(1)
      GOTO 120
  6   IF (FNUM(1) .LT. 0.0001) GOTO 120
      SCALE = FNUM(1)
      WRITE (LIS1, 125) SCALE
  125 FORMAT (' Scale from CONDA file: Scale =', F9.5)
      GOTO 120
  7   IF (NFNUM.LT.1) CALL KERNER (7, 'FCALIN')
      IF (FNUM(1).GT.0.0001) THEN
         BOV = FNUM(1)
         BP = BOV
         BR = BOV
         ENDIF
      IF (FNUM(2).GT.0.0001) BP = FNUM(2)
      IF (FNUM(3).GT.0.0001) BR = FNUM(3)
      WRITE (LIS1, 126) FNUM(1), FNUM(2), FNUM(3)
  126 FORMAT (' Temp. factors from CONDA: Bov=',
     *          F6.3, ' Bp =', F6.3, ' Br =', F6.3)
      IF (FNUM(2).LE.0.0001 .OR. FNUM(3).LE.0.0001)
     *   WRITE (LIS1, 127) BOV, BP, BR
  127 FORMAT (' Temp. factors to be used: Bov=',
     *          F6.3, ' Bp =', F6.3, ' Br =', F6.3)
      GOTO 120
  8   IF (NFNUM.NE.1) CALL KERNER (8, 'FCALIN')
      IF (FNUM(1) .LT. 0.0001 .OR. FNUM(1) .GE. STLMAX) GOTO 120
      STLMAX = FNUM(1)
      WRITE (LIS1, 128) STLMAX
  128 FORMAT (' Skip reflections if  sin(th)/lambda  >', F8.4)
      GOTO 120
  9   SMM =FNUM(1)
      IF(SMM.LT.0.0001) SMM=0.000001
      WRITE (LIS1, 129) SMM
  129 FORMAT (' Input factor (for testruns) SMM =     ', F8.4)
      GOTO 120
  10  J = 0
      DO 137 I = 1, 3
      IF (FNUM(I).LT.0.0001) GOTO 137
      K = NINT(FNUM(I))
      IF (K .GE. MAXHKL(I)) GOTO 137
      J = 1
      MAXHKL(I) = K
  137 CONTINUE
      IF (J .EQ. 0) GOTO 120
      WRITE (LIS1, 138) MAXHKL
  138 FORMAT (' Skip reflections if indices exceed MAXHKL =', 3I3)
      CALL KERI2F (MAXHKL, HKLMAX, 3)
      GOTO 120
  11  IOMAPS = 1
      GOTO 120
  12  IF (FNUM(1) .GE. 0.01) P77 = FNUM(1)
      GOTO 120
  140 CALL FILCLO (ICON, 'KEEP')
      SMAX = STLMAX
      NSET = 0
      CALL FILINQ (IATOMS, 'ATOMS', 'FORMATTED', 'INPUT', KINQ)
      IF (KINQ.EQ.-1) CALL KERROR (' No ATOMS file found',
     *   140, 'FCALIN')
      CALL KERINA (IATOMS, LIT, 1, LEND)
      IF (LIT(1) .NE. 'ATOMS') CALL KERROR
     *   (' Incorrect header on ATOMS file', 143, 'FCALIN')
      BACKSPACE IATOMS
      CALL ATOMIN (IATOMS, ATXYZ, ATNAME, IZAT, MAXAT, NAT, KEYT)
      NSET = NSET + 1
      IF (KEYT .GT. 1) WRITE (LIS1, FMT='(
     *  '' Temperature factors present in ATOMS file ignored'')')
      KEYT = 1
      NATQ = NAT
      NATH = 0
      N = 1
  143 CONTINUE
      IF (ATNAME(N)(1:1).EQ.'H' .AND. IZAT(N).EQ.1) NATH = NATH + 1
      IF (ATNAME(N)(1:1).EQ.'H' .OR.  IZAT(N).EQ.1 .OR.
     *      ATNAME(N)(1:1) .EQ. 'Q') THEN
         IF (N .EQ. NAT) GOTO 148
         DO 146 N1 = N, NAT - 1
         CALL KERNAB (ATXYZ(1,N1+1), ATXYZ(1,N1), 4)
         ATNAME(N1) = ATNAME(N1+1)
  146    IZAT(N1) = IZAT(N1+1)
  148    NAT = NAT - 1
         N = N - 1
         ENDIF
      CALL KERNZA (0.0, ATXYZ(5,N), 6)
      N = N + 1
      IF (N .LE. NAT) GOTO 143
      IF (NATH .NE. 0) WRITE (LIS1, FMT=
     *  '('' Number of H atoms rejected:'', I3)') NATH
      N = NATQ - NAT - NATH
      IF (N .GT. 0) WRITE (LIS1, FMT=
     *  '('' Nr of Q-atoms (= peaks) rejected:'', I3)') N
      IF (NAT .LE. 0) CALL KERROR ('.... No atoms left!', 0, 'FCALIN')
      CALL ATOMPR (LIS1, 7, ATXYZ, ATNAME, IZAT, NAT)
      CALL KERNAB (BUFFO(5), BUFFC(5), MAXBUF - 4)
      BUFFC(16) = SCALE
      BUFFC(17) = BOV
      BUFFC(18) = SCALE
      BUFFC(19) = BP
      BUFFC(20) = BR
      BUFFC(21) = STLMAX
      CALL KERNAB (HKLMAX, BUFFC(22), 3)
      CALL FCALCI (KEYT, ATXYZ, IZAT, ITAT, NAT)
      BUFFC(25) = NAT
      BUFFC(26) = P1SQ
      BUFFC(27) = PSQ
      CALL FCALCX
      RETURN
      END
      SUBROUTINE FCALCX
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      EQUIVALENCE (IPR1, IFILE(6)), (LIS1, IFILE(7)), (LIS2, IFILE(8))
      EQUIVALENCE (IBINFO, IFILE(11)), (IBINFC, IFILE(12))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     *               WAVE,     CELALL(10),  AMOLW,      ZET,
     *               NELEC,    F000,        ABSMU,      ICENT,
     *               ILATT,    ISYST,       ILAUE,      IMULT,
     *               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     *         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     *         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      COMMON /CRYSB/ SPGR,     WAVEAT,      CELATY(10)
      CHARACTER      SPGR *16, WAVEAT *2,   CELATY *2
      COMMON /FCALCA/ BP,       BR,       SCALE,    HKLMAX(3), STLMAX,
     *                IZTYPE(10), CELPAR(10), PSQ,  P1SQ,     ITRS(24),
     *        AMULT,  ASYMM,    ALATT,    ASYMCL,   NSYMC,    ASYMC,
     *                HKLX(3,24), IDHKL(24), HCODE, FOBS,     SIG,
     *                STL,      STL2,     ISS,      ENORM,
     *                FP,       PHIP,     FAP,      FBP,      EPSIL,
     *                EPSIL2,   SF2,      SF2P,     FPEXP(2,24)
      DIMENSION FITFO(3), FITFC2(51)
      EQUIVALENCE (HCODE, FITFO(1)), (EPSIL2, FITFC2(1))
      PARAMETER (MAXAT  = 993)
      PARAMETER (MAXBUF = 198)
      COMMON /  / SICO(12500), FF(500,10),  EXPBP(500),   EXPBR(500),
     *            SUMF2(500),  SUMF2P(500), SFAC(13,10),
     *            ATXYZ(10,MAXAT), IZAT(MAXAT), ITAT(MAXAT), NAT,
     *            BUFFO(MAXBUF),  BUFFC(MAXBUF), BUFBUF(MAXBUF)
      DIMENSION BLACOM(42000)
      EQUIVALENCE (BLACOM(1), SICO(1))
      COMMON /ATNAMA/ ATNAME(MAXAT)
      CHARACTER * 6   ATNAME
      WRITE (LIS1, FMT='('' Calculate structure factors'')')
      FICENT = FLOAT (ICENT)
      WRITE (LIS1, 183) P1SQ
  183 FORMAT (' Scattering fraction of known part:'/
     *   '    excluding symmetry related molecules: P1**2 =', F6.3)
      NITFC = 3 + 2 * NSYMM
      CALL BINOFF (27, IBINFC,'BINFC2', FITFC2, NITFC, BUFFC, KENDFC)
      NREFL  = 0
      IREFL  = 0
      SUMNR2 = 0.
      R2SUMA = 0.
      R2SUMB = 0.
      ALATT2 = NLATT * NLATT
      CALL BINIFF (1, IBINFO, 'BINFO', FITFO, NITFO, BUFFO, KENDFO)
      NITFO = 3
  200 CALL BINIFF (0, IBINFO, 'BINFO', FITFO, NITFO, BUFFO, KENDFO)
      IF (KENDFO.LT.0) GOTO 220
      IREFL = IREFL + 1
      CALL HKLC1U (HCODE, HKLX)
      CALL HKLSTL (HKLX, STL, STL2)
      IF (STL.GT.STLMAX .OR. HKLX(1,1).GT.HKLMAX(1) .OR.
     *      HKLX(2,1).GT.HKLMAX(2) .OR. HKLX(3,1).GT.HKLMAX(3)) THEN
         FITFC2(1) = -999.
         GOTO 215
         ENDIF
      CALL FCALX1 (ATXYZ, ITAT, NAT)
      SUMNR2 = SUMNR2 + ASYMM / EPSIL2
      SF2 = SUMF2 (ISS)
      SF2P= SUMF2P(ISS)
      NREFL = NREFL + 1
      G2 =  EXPBR(ISS) * EXPBR(ISS) * SUMF2(ISS) * ALATT2
      EO2 = ( SCALE * FOBS ) **2 / G2
      EO4 = EO2 * EO2
      DO 207 I = 1, NSYMM
      IF (FPEXP(1,I) .LT. 0.) GOTO 207
      R2SUMA = R2SUMA + (EO2 - FPEXP(1,I)*FPEXP(1,I)/G2)**2
      R2SUMB = R2SUMB + EO4
  207 CONTINUE
  215 CONTINUE
      CALL BINOFF (0, IBINFC, 'BINFC2',FITFC2, NITFC, BUFFC, KENDFC)
      GOTO 200
  220 CONTINUE
         CALL BINOFF (-1, IBINFC, 'BINFC2',FITFC2, NITFC, BUFFC, KENDFC)
      IF (IREFL.NE.NREFL) WRITE (LIS1, 242) NREFL
  242 FORMAT (' Note:', I6,' input reflections accepted ')
      WRITE (LIS1, 254) NINT(SUMNR2)
  254 FORMAT (' Number of reflections after expansion:', I12)
      WRITE (IPR1, 258)
  258 FORMAT (' Structure factor calculation finished')
      R2 = R2SUMA / R2SUMB
      R2E = 1. - P1SQ
      WRITE (LIS1, FMT='(/'' R2 for the model in SpGr P1:'', F5.2,
     *  ''   ( expected:'', F5.2, '' )''/ 34(''-'')/ )') R2, R2E
      RETURN
      END
      SUBROUTINE FCALX1  (ATXYZ, ITAT, NAT)
      DIMENSION ATXYZ(10,NAT), ITAT(NAT)
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      COMMON /FCALCA/ BP,       BR,       SCALE,    HKLMAX(3), STLMAX,
     *                IZTYPE(10), CELPAR(10), PSQ,  P1SQ,     ITRS(24),
     *        AMULT,  ASYMM,    ALATT,    ASYMCL,   NSYMC,    ASYMC,
     *                HKLX(3,24), IDHKL(24), HCODE, FOBS,     SIG,
     *                STL,      STL2,     ISS,      ENORM,
     *                FP,       PHIP,     FAP,      FBP,      EPSIL,
     *                EPSIL2,   SF2,      SF2P,     FPEXP(2,24)
      COMMON /  / SICO(12500), FF(500,10),  EXPBP(500),   EXPBR(500),
     *            SUMF2(500),  SUMF2P(500), SFAC(13,10)
      DIMENSION BLACOM(42000)
      EQUIVALENCE (BLACOM(1), SICO(1))
      DIMENSION FFF(10)
      S = STL * 400. + 1.
      IS = IFIX(S)
      STLDEL = S - FLOAT(IS)
      ISS = NINT(S)
      DO 110 J=1,NTYPE
      IF (CELPAR(J).LE.0.0) GOTO 110
      FFF(J) = FF(IS,J) + (FF(IS+1,J)-FF(IS,J)) * STLDEL
  110 CONTINUE
      CALL HKLEX1 (HKLX, HKLX)
      CALL HKLEX2 (HKLX, IDHKL, IEPS, IEPS2)
      EPSIL = IEPS
      EPSIL2 = IEPS2
      DO 600 J=1,NSYMM
      IF (IDHKL(J).EQ.0) GOTO 200
      K = IABS(IDHKL(J))
      FPEXP(1,J) = -K
      FPEXP(2,J) = FPEXP(2,K)
      IF (FPEXP(2,J).LT.0.0001) FPEXP(2,J)=0.0001
      IF (IDHKL(J).LT.0) FPEXP(2,J)=-FPEXP(2,J)
      GOTO 600
  200 FAP = 0.0
      FBP = 0.0
      DO 400 I=1,NAT
      TRIG = HKLX(1,J)*ATXYZ(1,I) + HKLX(2,J)*ATXYZ(2,I) +
     *       HKLX(3,J)*ATXYZ(3,I)
      IF (TRIG.LT.0.0) TRIG = TRIG - 0.00010
      ITRIG = MOD ( IFIX(TRIG * 10000. + 0.5), 10000)
      IF (ITRIG.LE.0) ITRIG = ITRIG + 10000
      IJ = ITAT(I)
      FAP = FAP + SICO(ITRIG + 2500) * FFF(IJ)
      FBP = FBP + SICO(ITRIG)        * FFF(IJ)
  400 CONTINUE
         FP = ALATT * SQRT (FAP*FAP + FBP*FBP) * EXPBP(ISS)
      PHIP = 0.0
      IF (FP.GT.0.001) PHIP = ATAN2(FBP,FAP) / 0.0174532925
      IF (PHIP.LT.0.0) PHIP = PHIP + 360.
      FPEXP(1,J) = FP
      FPEXP(2,J) = PHIP
  600 CONTINUE
      RETURN
      END
      SUBROUTINE TRACIN
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      EQUIVALENCE (ICON,   IFILE(4))
      EQUIVALENCE (LIS2,   IFILE(8))
      EQUIVALENCE (IBINFO, IFILE(11))
      EQUIVALENCE (IBINFC, IFILE(12))
      LOGICAL SWPRI
      EQUIVALENCE (SWPRI, SWITCH(10))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      COMMON /CRYSB/ SPGR,     WAVEAT,      CELATY(10)
      CHARACTER      SPGR *16, WAVEAT *2,   CELATY *2
      PARAMETER (KFA = 20671)
      PARAMETER (MAXBUF = 198)
      COMMON /  /  FA(KFA), FB(KFA),
     +             BUFFC2(MAXBUF), BUFFO(MAXBUF), BUFFFT(MAXBUF),
     +             FITFC2(51),     FITFO(3),      FITFFT(5),
     +             NITFC2, KENDFC, NITFO, KENDFO
      DIMENSION BLACOM(42000)
      EQUIVALENCE (BLACOM(1), FA(1))
      PARAMETER (MXP=50, MXP1=MXP + 1)
      COMMON /PIEK/ RHO(MXP1), SH(3,MXP1), IBOTS(MXP1), BOTS(MXP1), NPK
      COMMON /COTRA/ MULS(48,48), INDX(3,3), IDIMF, TT(3,48),
     *               EFST(48), FHST(48), PHST(48), CELLN(6), MAXHKL(3),
     *  P1SQ, PSQ, EMIN, EO2AV, SCX, BOVX, DAMPX, SMM, IOMAPS, SMAX,
     *  IRR(3,3,48), IORG, BBB, D2R, R2D,
     *  AMULT, ASYMM, ALATT, ASYMCL, NSYMC, ASYMC
      DIMENSION IRP(3,3)
      DIMENSION EO2A(4,2)
      DATA EO2A /1.000,1.250,2.000,3.250,1.000,1.5705,2.5251,3.9080/
      D2R = ATAN(1.0) / 45.0
      R2D = 1.0 / D2R
      CALL KERNZA(0.0, RHO, MXP1)
      CALL KERNZA(0.0, SH, 3*MXP1)
      NPK   = 0
      NR    = 0
      AMULT = FLOAT (IMULT)
      ASYMM = FLOAT (NSYMM)
      ALATT = FLOAT (NLATT)
      ASYMCL= FLOAT (ICENT*NLATT)
      NSYMC = NSYMM * ICENT
      ASYMC = FLOAT (NSYMC)
      CALL BINIFF (1,IBINFC,'BINFC2',FITFC2,NITFC2,BUFFC2,KENDFC)
      SC   = BUFFC2(18)
      BH   = BUFFC2(19)
      BOV  = BH
      PSQMAX = PSQ
      P1SQ = BUFFC2(26)
      PSQ  = BUFFC2(27)
      IF (PSQMAX .GE. 0.01) PSQ = AMIN1(PSQ, PSQMAX)
      SMAX = BUFFC2(21)
      CALL KERNAB (BUFFC2, BUFFFT, 27)
      BUFFFT(28) = 5.
      CALL BINIFF (1,IBINFO,'BINFO',FITFO,NITFO,BUFFO,KENDFO)
      IEMIN = IFIX (EMIN * 2. + 1.5)
      IF (IEMIN.GT.4) IEMIN=4
      EMINX = FLOAT(IEMIN)/2. -0.5
      IF (ABS(EMINX - EMIN) .GT. 0.01) THEN
         EMIN = EMINX
      WRITE (CHOUT, FMT = '('' Reset value for Emin: '', F10.2)') EMIN
         CALL SHOUT
         ENDIF
      EO2AV = EO2A(IEMIN,ICENT)
      BBB = EO2AV * (1. - PSQ * P1SQ) - (PSQ - PSQ * P1SQ)
      IF (SWPRI) WRITE (LIS2, 334) P1SQ, PSQ, EMIN, EO2AV, BBB
  334 FORMAT (4X, 'P1SQ       PSQ       EMIN   ',
     +       '   EO2AV      BBB' /  5F10.5  )
      IF (SWPRI) WRITE (LIS2, 336) SCX, BOVX, DAMPX, SMM
  336 FORMAT (3X,'Scale',6X,'BH',6X,'  Damp-F',6X,
     *      'SMM  sharpening params' / 6F10.5 )
      DO 402 I=1,NSYMM
      L=I+NSYMM
      DO 402 J=1,3
      TT(J,I)=TSYMM(J,I)
      IF(ICENT.EQ.2) THEN
         TTJL=AMOD(1.0-TT(J,I),1.0)
         TT(J,L)=TTJL
         ENDIF
      DO 402 K=1,3
      IF(ICENT.EQ.2) IRR(K,J,L) = -IRSYMM(K,J,I)
  402 IRR(K,J,I) = IRSYMM(K,J,I)
      WRITE (LIS2, FMT='('' Symmetry operations for sp.gr. '', A16)')
     *   SPGR
      DO 419 I=1,NSYMC
      WRITE (LIS2, 407) I, ((IRR(L,M,I),M=1,3), TT(L,I), L=1,3)
  407 FORMAT (I4, ' :', 3('  / ', 3I3, F8.3))
      DO 419 J=1,NSYMC
      CALL MATAXI (IRR(1,1,I), IRR(1,1,J), IRP)
      DO 418 M=1,NSYMC
      IFIT=0
      DO 415 K=1,3
      DO 415 L=1,3
      IF (IRR(K,L,M) .EQ. IRP(K,L)) IFIT=IFIT+1
  415 CONTINUE
      IF (IFIT .NE. 9) GOTO 418
      MULS(I,J)=M
      GOTO 419
  418 CONTINUE
      CALL KERROR ('Symmetry error' , 0, 'SYTRA')
  419 CONTINUE
      IF (SWPRI)  THEN
         WRITE (LIS2,722) (I,I=1,NSYMC)
  722    FORMAT('0Multiplication table for symmetry operations'/
     +       (5X,24I3/))
         DO 730 I=1,NSYMC
  730    WRITE (LIS2, 732) I, (MULS(I,J),J=1,NSYMC)
  732    FORMAT(I3, 2X, 24I3)
         ENDIF
      RETURN
      END
      SUBROUTINE COTRA1
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      EQUIVALENCE (IPR1, IFILE(6)), (LIS1, IFILE(7)), (LIS2, IFILE(8))
      EQUIVALENCE (IBINFO, IFILE(11)), (IBINFC, IFILE(12))
      EQUIVALENCE (ITS14,  IFILE(14)), (ITS15,  IFILE(15))
      EQUIVALENCE (IBINFF, IFILE(16))
      LOGICAL SWPRI
      EQUIVALENCE (SWPRI, SWITCH(10))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      PARAMETER (KFA = 20671)
      PARAMETER (MAXBUF = 198)
      COMMON /  /  FA(KFA), FB(KFA),
     +             BUFFC2(MAXBUF), BUFFO(MAXBUF), BUFFFT(MAXBUF),
     +             FITFC2(51),     FITFO(3),      FITFFT(5),
     +             NITFC2, KENDFC, NITFO, KENDFO
      DIMENSION BLACOM(42000)
      EQUIVALENCE (BLACOM(1), FA(1))
      DIMENSION  BUF(2,24)
      EQUIVALENCE (EPSIL2, FITFC2(1)), (SFTOT2,  FITFC2(2)),
     +            (SFPAR2, FITFC2(3)), (BUF(1,1),FITFC2(4))
      PARAMETER (MXP=50, MXP1=MXP + 1)
      COMMON /PIEK/ RHO(MXP1), X1(3,MXP1), IBOTS(MXP1), BOTS(MXP1), NPK
      COMMON /COTRA/ MULS(48,48), INDX(3,3), IDIMF, TT(3,48),
     *               EFST(48), FHST(48), PHST(48), CELLN(6), MAXHKL(3),
     *  P1SQ, PSQ, EMIN, EO2AV, SCX, BOVX, DAMPX, SMM, IOMAPS, SMAX,
     *  IRR(3,3,48), IORG, BBB, D2R, R2D,
     *  AMULT, ASYMM, ALATT, ASYMCL, NSYMC, ASYMC
      DIMENSION SHKL(3,24), ISHKL(3,48), LHK(3), M3D(3), IORGX(14)
      DIMENSION MHKL(3), M2HKL(3), KM3D(3)
      LOGICAL SWHKL, SWHKL2
      DIMENSION MSS(4)
      LOGICAL  LMSS(4), LMSS4, LMSS5, LMSS6, LMSS7
      EQUIVALENCE (LMSS(1), LMSS4), (LMSS(2), LMSS5),
     *            (LMSS(3), LMSS6), (LMSS(4), LMSS7)
      EQUIVALENCE (MSS(1), MSS4), (MSS(2), MSS5),
     *            (MSS(3), MSS6), (MSS(4), MSS7)
      DATA ZERO  /0.0 / , IZRO, NNN / 0, 999/
      DATA IORGX /1,1,1,2,    2, 5,  5,  3,  3,  3,    3,    3,  4,  4/
      DATA I1, I2 / 0, 0/
      IORG = IORGX(ILAUE)
      CALL KERNZI (0, INDX, 9)
      CALL KERNZI (0, KM3D, 3)
      CALL KERNZI (0, M2HKL, 3)
      CALL KERNAI (MAXHKL, MHKL, 3)
      CALL KERNZI (0, MAXHKL, 3)
      IF (IPOLA .NE. 8) GOTO 102
      DO 101 I = 1, 3
      CELLN(I) = 1.4 * CELL(I)
      CELLN(I+3) = 90.
      M2HKL(I) = 2 * MHKL(I)
  101 M3D(I) = 6. * SMAX * CELL(I)
      CELLN(6) = 120.
      GOTO 110
  102 DO 103 I=1,3
      M3D(I) = SMAX * 2.31 * CELL(I)
      M2HKL(I) = MHKL(I)
      IF (ILAUE .EQ. 3) M3D(I) = SMAX * 2.0 * CELL(I)
      IF (IUNIQ .EQ. I) M3D(I) = SMAX * 2.0 * CELL(I)
      CELLN(I+3)= CELL(I+3)
  103 CELLN(I)  = CELL(I)/2.0
      GOTO (110,104,108,104,104), IORG
  104 CELLN(1)  = CELL(1)
      CELLN(2)  = CELLN(1)
      M3D(1) = SMAX * 4.0 * CELL(1)
      M3D(2) = M3D(1)
      M3DMMM = SQRT ( 0.5 * FLOAT( KFA - 300) )
      KM3D(1) = MIN0 (M3D(1), M3DMMM)
      KM3D(2) = KM3D(1)
      KM3D(3) = MHKL(3)
      M2HKL(1) = 2 * MHKL(1)
      M2HKL(2) = 2 * MHKL(1)
      IF (IORG .EQ. 2) GOTO 110
      CELLN(3)  = CELL(1)
      M3D(3) = M3D(1)
      M2HKL(3) = MHKL(1)
      IF (IORG .EQ. 4) GOTO 110
      M3D(1) = SMAX * 4.62 * CELL(1)
      M3D(2) = M3D(1)
      M3D(3) = M3D(1)
      M2HKL(1) = 2 * MHKL(1)
      M2HKL(2) = 2 * MHKL(1)
      M2HKL(3) = 2 * MHKL(1)
      GOTO 110
  108 CELLN(1)  = CELL(1)
      CELLN(2)  = CELLN(1)
      M3D(1) = SMAX * 4.62 * CELL(1)
      M3DMMM = SQRT ( 0.5 * FLOAT( KFA - 300) )
      M3D(1) = MIN0 (M3D(1), M3DMMM)
      M3D(2) = M3D(1)
      M2HKL(1) = 2 * MHKL(1)
      M2HKL(2) = 2 * MHKL(1)
      M2HKL(3) = 2 * MHKL(1)
  110 WRITE (LIS2, 112) CELL, VOLUM, CELLN, IORG, M3D, M2HKL, KM3D
  112 FORMAT (/' CELL: ', 3F8.3, 3F7.2, ' VOLUME:' ,F9.2/
     *         ' CELLN:', 3F8.3, 3F7.2, ' transformed cell'/
     *         ' (Origin keys: IORG =',I2, ')' /
     *         ' (Expected maximum transformed indices:  ', 3I3,' )'/
     *         ' ( [ TEMP: M2HKL=',3I3,'    TEMP: KM3D=', 3I3,' )' /)
      IF (IPOLA .EQ. 0) THEN
         IDIMF=3
         CALL BINOFF (28, IBINFF, 'BINFFT', FITFFT, 5, BUFFFT, KENDFF)
         GOTO 130
         ENDIF
      WRITE (LIS2, 119)
  119 FORMAT(' POLAR SPACE GROUP - NO 3-D SEARCH')
      IDIMF=2
      IF (IPOLA.EQ.3 .OR. IPOLA.EQ.5 .OR. IPOLA.EQ.6) IDIMF=1
      CALL FILINQ (ITS15, 'SCRA15', 'UNFORMATTED', 'SCRATCH', KINQ)
      GOTO (122, 120, 123, 121, 122, 121, 120, 121), IPOLA
  120 INDX(1,1) = 1
      INDX(2,3) = 1
      INDX(3,2) =-1
      M3D(2)=0
      GOTO 130
  121 I=1
      GOTO 125
  122 I=2
      GOTO 125
  123 I=3
  125 J=MOD(I,3)+1
      K=6-I-J
      INDX(1,I) = 1
      INDX(2,J) = 1
      INDX(3,K) = 1
      M3D(K)=0
      IF(IDIMF.EQ.1) M3D(J)=0
  130 CONTINUE
      CALL KERNZA (0.0, FA, KFA)
      CALL KERNZA (0.0, FB, KFA)
      KLAAR = 0
  131 M2 = 2 * M3D(3) + 1
      M1 = M2 * (2 * M3D(2) + 1)
      M3 = M3D(2) * M2
      NEED = M3D(1) * M1 + M3D(2) * M2 + M3D(3)
      WRITE (LIS2, 132) M3D,NEED
  132 FORMAT(' Max indices and storage for general scan: ',3I3,I8)
      M5 = 2 * M3D(2) + 1
      MS1 = (M3D(2) + 1) * M2
      MS2 = (M3D(1) + 1) * M2
      MS3 = (M3D(1) + 1) * M5
      MS12= MS1 + MS2
      MS123= MS12 + MS3
      KFA1 = KFA - 100
      IF (MS123 .GT. KFA1)
     *   WRITE (LIS2, 133)  MS1, MS2, MS3, MS123, KFA1
  133    FORMAT(/' Storage 2D scans:    BC    AC    AB    Sum   LIMIT:'/
     *           ' -----------------', 3I6, I7, I8 /)
      KLAAR = KLAAR + 1
      IF (MS123 .GT. KFA1 .AND. KLAAR .LE. 4) THEN
         KLA = NINT ( 100. * SQRT ( FLOAT( MS123) / FLOAT( KFA1 ) ) )
         DO 134 I = 1, 3
         M3D(I) = 100 * M3D(I) / KLA
  134    CONTINUE
         GOTO 131
         ENDIF
      IF (MS123 .GT. KFA1) THEN
         WRITE (LIS1, 133)  MS1, MS2, MS3, MS123, KFA1
         WRITE (LIS1, FMT='(/'' Rerun TRACOR with lower sinth/lam'')')
         CALL KERROR ('Storage problems: KFA1.', 133, 'COTRA1')
         ENDIF
      MSTOT = MS123
      NKFA = MIN0 (NEED, KFA)
      IDIMFX = IDIMF
      IF (ILAUE .EQ. 1) IDIMFX = 4
      IF (NEED.LE.KFA .AND. IDIMFX.EQ.3) IDIMFX = 5
      MSS4 = MS123 + MS2
      MSS5 = MSS4  + MS2
      MSS6 = MSS5  + MS2
      MSS7 = MSS6  + MS2
      DO 135 M = 1,4
      LMSS(M) = .TRUE.
      IF (MSS(M) .LE. KFA - 100) THEN
         LMSS(M) = .FALSE.
         MSTOT = MSS(M)
         ENDIF
  135 CONTINUE
      NFOUR = 0
      KBOT = 1
      KTOP = NEED
      NPASS = 1
      IPASS = 0
      IF (NEED .LE. KFA .OR. ILAUE .EQ. 1) GOTO 140
      CALL FILINQ (ITS14, 'SCRA14', 'UNFORMATTED', 'SCRATCH', KINQ)
      KTOP = KFA - MSTOT
      NPASS = (NEED - 1 + MSTOT) / KFA + 1
  140 NACC  = 0
      NRIN  = 0
      NRINP = 0
      SEOBS = 0.0
      SES2  = 0.0
      TOTREF= 0.0
      NUNOBS = 0
      IADDM = -9999
      NREF14  =  0
      SWHKL = .FALSE.
      SWHKL2 = .FALSE.
      IF (.NOT. SWPRI) SWHKL2 = .TRUE.
      ITREF = 0
      ITREFE = 0
      ITREFS = 0
      IF (SWPRI) WRITE (LIS2, 150)
  150   FORMAT('0  h   k   l   Eobs  eps   SF2   SF2P',
     *  '   Fp   phase   Fp   phase   .....')
  160 CALL BINIFF (0,IBINFC,'BINFC2',FITFC2,NITFC2,BUFFC2,KENDFC)
      IF (KENDFC .LT. 0) GOTO 500
      NSREF = NSYMM / NINT (EPSIL2)
      NRIN = NRIN + NSREF
      NRINP = NRINP + 1
      CALL BINIFF (0,IBINFO,'BINFO',FITFO,NITFO,BUFFO,KENDFO)
      IF (KENDFO .LT. 0) CALL KERNER(-2, 'COTRA1')
      IF (FITFC2(1) .LT. 0.) THEN
         ITREFS = ITREFS + 1
         GOTO 160
         ENDIF
      CALL HKLC1U (FITFO(1), SHKL)
      CALL HKLSTL (SHKL, STL, STL2)
      IF (STL .GT. SMAX) GOTO 160
      IF (FITFO(2) .LE. 6.01*FITFO(3)) THEN
         NUNOBS = NUNOBS + 1
         FITFO(2) = FITFO(2) / 4.
         ENDIF
      EF = SCX * FITFO(2) * EXP(BOVX * STL2) / SQRT(SFTOT2 * ALATT)
      IF (EF.GT.3.0) EF=3.0
      IF (EF.LT.EMIN) THEN
            ITREFE = ITREFE + 1
            GOTO 160
            ENDIF
      NACC = NACC + NSREF
      SXX=STL/SMM
      IF(SXX.GT.SMM) SXX=1.0
      SXX = SQRT(SXX) * EXP(DAMPX * STL2) / SQRT(SFPAR2 * ALATT)
      CALL HKLEX1 (SHKL,  SHKL)
      CALL HKLEX2 (SHKL, ISHKL, IEPS, IEPS2)
      EPSIL = IEPS
      CALL KERF2I (SHKL, ISHKL, 3*NSYMM)
      IF (SWHKL2) GOTO 172
      IF (SWHKL) THEN
         SWHKL2 = .TRUE.
         SWHKL = .FALSE.
         GOTO 172
         ENDIF
      IF (ISHKL(1,1).EQ.0 .OR. ISHKL(2,1).EQ.0 .OR.
     *    ISHKL(3,1).EQ.0) GOTO 172
      IF (ISHKL(1,1).EQ.ISHKL(2,1) .OR. ISHKL(2,1).EQ.ISHKL(3,1) .OR.
     *    ISHKL(3,1).EQ.ISHKL(1,1)) GOTO 172
      IF (ISHKL(1,1).EQ.-ISHKL(2,1) .OR. ISHKL(2,1).EQ.-ISHKL(3,1) .OR.
     *    ISHKL(3,1).EQ.-ISHKL(1,1)) GOTO 172
      SWHKL = .TRUE.
  172 IF (SWPRI .AND. (SWHKL .OR. NRINP.LT.22)) WRITE
     *   (LIS2,180) (ISHKL(IJJ,1),IJJ=1,3), EF, (FITFC2(I),I=1,NITFC2)
  180 FORMAT (3I4, F5.2, F3.0, 2F7.0, 1X, 7(F7.2,F6.1)
     *        / 35X, 7(F7.2,F6.1))
      SUMES=0.0
      DO 190 IS=1,NSYMM
      PHST(IS) = BUF(2,IS)
      IF (PHST(IS) .LT. 0.0) PHST(IS) = 360. + PHST(IS)
      IF (BUF(1,IS) .GE. 0.0) THEN
         EFST(IS) = BUF(1,IS) * SXX
         SUMES = SUMES + EFST(IS)**2
      ELSE
         ISK = NINT (-BUF(1,IS))
         EFST(IS) = EFST(ISK)
         ENDIF
  190 CONTINUE
      IF (ICENT.EQ.1) GOTO 196
      DO 195 IS = NSYMM+1 , NSYMC
      DO 194 I=1,3
  194 ISHKL(I,IS)= - ISHKL(I,IS-NSYMM)
      EFST(IS) = EFST(IS-NSYMM)
  195 PHST(IS) = 360. - PHST(IS-NSYMM)
  196 CONTINUE
      TOTREF = TOTREF + 1.
      SUMES = SUMES / FLOAT(NSREF)
      SES2 = SES2 + SUMES
      EF2 = EF**2
      SEOBS = SEOBS + EF2
      EF2 = EF2 - PSQ * SUMES - BBB
      DO 400 I=2,NSYMC
      DO 390 II=1,NSYMM
      PTRA = SHKL(1,II)*TT(1,I)+ SHKL(2,II)*TT(2,I)+ SHKL(3,II)*TT(3,I)
      IJ = MULS(II,I)
      PTRA = PHST(II) - PHST(IJ) - 360.*PTRA
      EFTRA = EF2 * EFST(II) * EFST(IJ) / EPSIL
      IF (EFTRA .LT. 0.) THEN
         EFTRA = ABS(EFTRA)
         PTRA = PTRA + 180.
         ENDIF
      PTRA = AMOD(-PTRA, 360.)
      DO 301 J=1,3
  301 LHK(J) = ISHKL(J,II) - ISHKL(J,IJ)
      GOTO (302, 303, 303, 305, 304), IORG
  302 LHK(1)=LHK(1)/2
      LHK(2)=LHK(2)/2
  303 LHK(3)=LHK(3)/2
      GOTO 305
  304 IF (IPOLA .NE. 8) GOTO 305
      LHK(1) = LHK(1) - LHK(2)
      LHK(2) = LHK(2) - LHK(3)
      LHK(3) = 0
  305 DO 306 J=1,3
      IF (IABS(LHK(J)) .LE. M3D(J)) GOTO 306
      ITREF = ITREF + 1
      GOTO 390
  306 CONTINUE
      DO 308 J=1,3
  308 MAXHKL(J) = MAX0 (MAXHKL(J), IABS(LHK(J)))
      IF(LHK(1)) 320, 310, 330
  310 IF(LHK(2)) 320, 315, 330
  315 IF(LHK(3)) 320, 390, 330
  320 DO 322 J=1,3
  322 LHK(J)=-LHK(J)
      PTRA = - PTRA
  330 IADD = LHK(3) + LHK(2) * M2 + LHK(1) * M1
      IF (PTRA .LT. 0.) PTRA = PTRA + 360.
      IF (SWHKL) WRITE (LIS2, 331) (ISHKL(L,I),L=1,3),
     *   (ISHKL(L,II),L=1,3), (ISHKL(L,IJ),L=1,3), LHK, IADD
  331 FORMAT (' HKL for: I,II,IJ,M; IADD: ', 4(3I3,2X), I6)
      GOTO (338, 338, 343, 341, 338), IDIMFX
      CALL KERROR ('Impossible ...', 338, 'COTRA1')
  338 FA(IADD) = FA(IADD) + EFTRA * COS(PTRA * D2R)
      FB(IADD) = FB(IADD) + EFTRA * SIN(PTRA * D2R)
      GOTO 390
  341 FITFFT(1) = LHK(1)
      FITFFT(2) = LHK(2)
      FITFFT(3) = LHK(3)
      FITFFT(4) = EFTRA
      FITFFT(5) = PTRA
      CALL BINOFF(0, IBINFF, 'BINFFT', FITFFT, 5, BUFFFT, KENDFF)
      GOTO 390
  343 CONTINUE
      IF (LHK(1) .EQ. 0) THEN
         IADD  = LHK(2) * M2 + LHK(3)
         GOTO 338
         ENDIF
      IF (LHK(2) .EQ. 0) THEN
         IADD = MS1 + LHK(1) * M2 + LHK(3)
         GOTO 338
         ENDIF
      IF (LHK(3) .EQ. 0) THEN
         IADD = MS12 + LHK(1) * M5 + LHK(2)
         GOTO 338
         ENDIF
      IF (ILAUE .LE. 3) GOTO 341
      IF (LMSS4) GOTO 345
      IF (LHK(1) .EQ. LHK(2)) THEN
         IADD = MS123 + LHK(1) * M2 + LHK(3)
         GOTO 338
         ENDIF
      IF (LMSS5) GOTO 345
      IF (LHK(1) .EQ. -LHK(2)) THEN
         IADD = MSS4  + LHK(1) * M2 + LHK(3)
         GOTO 338
         ENDIF
      IF (LMSS6) GOTO 345
      IF (LHK(2) .EQ. -2 * LHK(1)) THEN
         IADD = MSS5  + LHK(1) * M2 + LHK(3)
         GOTO 338
         ENDIF
      IF (LMSS7) GOTO 345
      IF (2 * LHK(2) .EQ. -LHK(1)) THEN
         IADD = MSS6  - LHK(2) * M2 + LHK(3)
         GOTO 338
         ENDIF
  345 IF (IADD .GT. KTOP) GOTO 380
      IADD = IADD + MSTOT
      GOTO 338
  380 WRITE (ITS14) IADD, EFTRA, PTRA
      IADDM = MAX0(IADDM,IADD)
      NREF14 = NREF14 + 1
  390 CONTINUE
  400 CONTINUE
      GOTO 160
  500 IF (ITREF .GT. 0) WRITE (LIS1, FMT='(
     *   '' Nr of refl. bypassed (HKL > M3D)'', I6)') ITREF
      IF (ITREFE .GT. 0) WRITE (LIS1, FMT='(
     *   '' Nr of refl. bypassed (E < EMIN) '', I6)') ITREFE
      IF (ITREFS .GT. 0) WRITE (LIS1, FMT='(
     *   '' Nr of refl. bypassed (S > SMAX) '', I6)') ITREFS
      IF (NREF14 .GT. 0) THEN
         WRITE(LIS2, FMT = '('' max. indices-address      : '', I10,
     *        '' (calc:'',I10,'')'',/,
     *        '' Nr. of refl. to file ITS14: '', I10,/)')
     *        IADDM, NEED, NREF14
         NEED = IADDM
         WRITE (LIS1, FMT='('' Max. transformed hkl: '', 3I3)') MAXHKL
         ENDIF
      WRITE(LIS2, FMT = '('' Nr. of weak reflections  : '', I10)')
     *   NUNOBS
      WRITE (LIS2, FMT='('' Max. transformed hkl: '', 3I3)') MAXHKL
      IF (NPASS .GT. 1) WRITE (IPR1, FMT= '('' Approximately'', I3,
     *   '' passes needed for output/input scratch file: '')') NPASS
      ITREF = 0
      GOTO 630
  600 IDIMFX = 6
      KBOT=KTOP+1
      KTOP = MIN0 (KTOP + KFA, NEED)
      NKFA = KTOP - KBOT + 1
      MSTOT = 0
      CALL KERNZA (0.0, FA, NKFA)
      CALL KERNZA (0.0, FB, NKFA)
      REWIND ITS14
      NREFFI = 0
      NREFFO = 0
  620 READ (ITS14) IADD, F, P
      NREFFI = NREFFI + 1
      IF (IADD.LT.1) GOTO 625
      IF (IADD.LT.KBOT .OR. IADD.GT.KTOP) GOTO 620
      NREFFO = NREFFO + 1
      J = IADD - KBOT + 1
      FA(J) = FA(J) + F * COS(P * D2R)
      FB(J) = FB(J) + F * SIN(P * D2R)
      GOTO 620
  625 CONTINUE
      IF (IPASS .EQ. 0) WRITE(LIS2, FMT =
     *   '('' Nr. of refl. from file ITS14:'',I9 )') NREFFI
      WRITE(LIS2, FMT =
     *   '('' Nr. of reflections stored   :'',I9 )') NREFFO
  630 IF (IDIMFX .EQ. 4) GOTO 650
      DO 640 K = 1, NKFA
      J = K + KBOT - 1
      F = SQRT (FA(K)**2 + FB(K)**2)
      IF (F .LT. 0.000001) GOTO 640
      IF (F .LT. 0.01) THEN
         ITREF = ITREF + 1
         GOTO 640
         ENDIF
      NFOUR = NFOUR + 1
      GOTO (632, 632, 635,  631, 632, 632), IDIMFX
  631 CALL KERROR ('Impossible ...', 631, 'COTRA1')
  632 IADD = J
  633 I1 = (IADD + M3 + M3D(3)) / M1
      IADD  = IADD - I1 * M1
      I2 = (IADD + M3 + M3D(3)) / M2 - M3D(2)
      I3 = IADD - I2 * M2
      IF (IDIMF.EQ.3) GOTO 639
      I = I1 * INDX(1,1) + I2 * INDX(1,2) + I3 * INDX(1,3)
      I2= I1 * INDX(2,1) + I2 * INDX(2,2) + I3 * INDX(2,3)
      I1=I
      IF (IDIMF.EQ.1) WRITE(ITS15) I1,FA(K),FB(K)
      IF (IDIMF.EQ.2) WRITE(ITS15) I1,I2,FA(K),FB(K)
      GOTO 640
  635 IF (J. GT. MSTOT) GOTO 638
      IF (J .LE. MS1) GOTO 632
      IF (J .LE. MS12) THEN
         I1 = (J - MS1 + M3D(3)) / M2
         I2 = 0
         I3 = J - MS1 - I1 * M2
         GOTO 639
         ENDIF
      IF (J .LE. MS123) THEN
         I1 = (J - MS12 + M3D(2)) / M5
         I2 = J - MS12 - I1 * M5
         I3 = 0
         GOTO 639
         ENDIF
      IF (ILAUE .LE. 3) CALL KERROR (' Kan niet ', 635, 'COTRA')
      IF (J .LE. MSS4) THEN
         I1 = (J - MS123 + M3D(3)) / M2
         I2 = I1
         I3 =  J - MS123 - I1 * M2
         GOTO 639
         ENDIF
      IF (J .LE. MSS5) THEN
         I1 = (J - MSS4 + M3D(3)) / M2
         I2 = -I1
         I3 = J - MSS4 - I1 * M2
         GOTO 639
         ENDIF
      IF (J .LE. MSS6) THEN
         I1 = (J - MSS5 + M3D(3)) / M2
         I2 = -2 * I1
         I3 = J - MSS5 - I1 * M2
         GOTO 639
         ENDIF
      IF (J .LE. MSS7) THEN
         I2 = (J - MSS6 + M3D(3)) / M2
         I1 = 2 * I2
         I3 = J - MSS6 - I2 * M2
         I2 = - I2
         GOTO 639
         ENDIF
      CALL KERROR (' Kan niet ', 635, 'COTRA')
  638 IADD = J - MSTOT
      GOTO 633
  639 FITFFT(1) = I1
      FITFFT(2) = I2
      FITFFT(3) = I3
      FITFFT(4) = F
      P = ATAN2(FB(K),FA(K)) * R2D
      IF (P .LT. 0.0) P = P + 360.
      FITFFT(5) = P
      CALL BINOFF(0, IBINFF, 'BINFFT', FITFFT, 5, BUFFFT, KENDFF)
  640 CONTINUE
      IF (NREF14 .EQ. 0) GOTO 650
      IPASS = IPASS + 1
      IF (NPASS .GT. 1) THEN
         IF (IPASS .EQ. 1) WRITE (ITS14) IZRO,ZERO,ZERO
         WRITE (IPR1, FMT='('' Pass'', I3, '' completed'')') IPASS
         WRITE (LIS2, FMT='('' Pass'', I3, '' completed'')') IPASS
         IF (KTOP .LT. NEED) GOTO 600
         ENDIF
  650 IF (IDIMF.EQ.1) WRITE(ITS15) NNN,ZERO,ZERO
      IF (IDIMF.EQ.2) WRITE(ITS15) NNN,IZRO,ZERO,ZERO
      IF (IDIMF.NE.3) THEN
         REWIND ITS15
      ELSE
         CALL BINOFF(-1, IBINFF, 'BINFFT', FITFFT, 5, BUFFFT, KENDFF)
         ENDIF
      IF (ITREF .GT. 0) WRITE (LIS1, FMT='(
     *   '' Nr of refl. bypassed (F < 0.01) '', I6)') ITREF
      IF (SWPRI) THEN
         WRITE (LIS2, 902) NRIN,NACC, NINT(TOTREF), NFOUR
  902    FORMAT (' Nr of reflections input:', I9 /
     *   ' Nr of reflections accepted:', I6,' with', I6, ' independent.'
     *  /' Nr of terms that contributed to the Fourier map:', I8)
         WRITE (LIS2, 903) MAXHKL
  903    FORMAT (' Maximum indices for Fourier synthesis: ', 3I3)
         ENDIF
      CALL FILCLO (ITS14, 'DELETE')
      SEOBS = SEOBS / TOTREF
      SES2 = SES2 / TOTREF
      PPES = P1SQ * EO2AV + (1.-P1SQ)
      EOOR = EO2AV - PSQ * PPES - BBB
      ENUL = SEOBS - PSQ * SES2 - BBB
      WRITE (LIS2, 910) EO2AV,SEOBS,PPES,SES2,EOOR,ENUL
  910 FORMAT (/' Average of        theory       observed' /
     +' EOBS**2       ',F10.4,5X,F10.4 /
     +' ES**2         ',F10.4,5X,F10.4 /
     +' EO2-IMV-OR    ',F10.4,5X,F10.4 //)
      RETURN
      END
      SUBROUTINE D12FOU
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      EQUIVALENCE (LIS1,  IFILE(7))
      EQUIVALENCE (ITS15, IFILE(15))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      PARAMETER (NPKM2=100,MNGR2=163,         MNGR1=200,LC24=24)
      COMMON / / SICO(1250),   RHOSUM(MNGR2,MNGR2),    RHOMAX(NPKM2),
     +            IC1(MNGR1), IC2(MNGR1),  IMAX(NPKM2),   JMAX(NPKM2),
     +            XMAX(NPKM2), YMAX(NPKM2),  RHOP(NPKM2),  IRHO(LC24)
      DIMENSION BLACOM(42000)
      EQUIVALENCE (BLACOM(1), SICO(1))
      PARAMETER (MXP=50, MXP1=MXP + 1)
      COMMON /PIEK/ RHO(MXP1), SH(3,MXP1), IBOTS(MXP1), BOTS(MXP1), NPK
      COMMON /COTRA/ MULS(48,48), INDX(3,3), IDIMF, TT(3,48),
     *               EFST(48), FHST(48), PHST(48), CELLN(6), MAXHKL(3),
     *  P1SQ, PSQ, EMIN, EO2AV, SCX, BOVX, DAMPX, SMM, IOMAPS, SMAX,
     *  IRR(3,3,48), IORG, BBB, D2R, R2D,
     *  AMULT, ASYMM, ALATT, ASYMCL, NSYMC, ASYMC
      DIMENSION AL(3)
      CALL SICOT (SICO, 1250)
      IF (IPOLA .EQ. 8) THEN
         CALL KERNAB (CELLN, AL, 3)
         ANGLE = 120.
         GOTO 135
         ENDIF
      DO 112 L=1,2
      AL(L)=0.
      DO 110 I=1,3
      X=0.
      DO 108 J=1,3
  108 X = X + INDX(L,J) * RRMAT(I,J)
  110 AL(L) = AL(L) + X * INDX(L,I)
  112 AL(L)=SQRT(AL(L))
      AL(3)=0.
      DO 122 I=1,3
      X=0.
      DO 121 J=1,3
  121 X = X + (INDX(1,J) - INDX(2,J)) * RRMAT(I,J)
  122 AL(3) = AL(3) + X * (INDX(1,I) - INDX(2,I))
      AL(3)=(AL(1)*AL(1)+AL(2)*AL(2)-AL(3))/AL(1)/AL(2)
      ANGLE = ACOS(0.5 * AL(3)) * R2D
      IF (IDIMF .NE. 1) GOTO 135
      NGR1=AL(1)/0.21
      IF (NGR1.GT.200) NGR1=200
      WRITE (LIS1, 127) AL(1),NGR1
  127 FORMAT('0 FOURIER IS CALCULATED ALONG AXIS OF',F8.3,' A' /
     +' NUMBER OF GRID POINTS IN THIS DIRECTION (QX):',I5 )
      CALL D1FOUR (IOMAPS,NGR1,SCALE)
      GOTO 137
  135 NGR1=AL(1)/0.25
      IF (NGR1.GT.163) NGR1=163
      NGR2=AL(2)/0.25
      IF (NGR2.GT.163) NGR2=163
      WRITE (LIS1, 136) AL(1), AL(2), ANGLE, NGR1, NGR2
  136 FORMAT('0FOURIER IS CALCULATED IN PLANE OF AXES OF', /, 6X,
     +F8.3,' AND',F8.3,' A, ANGLE = ',F8.3 /
     +' NUMBER OF GRID POINTS IN FIRST  DIRECTION (QX):',I5/
     +' NUMBER OF GRID POINTS IN SECOND DIRECTION (QY):',I5)
      CALL D2FOUR (IOMAPS,NGR1,NGR2,SCALE)
  137 WRITE (LIS1, 138) SCALE
  138 FORMAT(' FOURIER SCALEFACTOR = ',F10.4)
      CALL FILCLO (ITS15, 'DELETE')
      CALL KERNZA (0.0, SH, 3 * MXP1)
      RHOMIN = 0.5 * RHOP(1)
      DO 206 I=1,NPK
      RHO(I) = RHOP(I)
      IF (RHO(I) .LT. RHOMIN) THEN
         NPK=I-1
         RETURN
         ENDIF
      DO 204 K=1,3
      IF (INDX(1,K) .NE. 0) SH(K,I) = XMAX(I)
  204 CONTINUE
      IF (IDIMF .EQ. 1) GOTO 206
      DO 205 K=1,3
      IF (INDX(2,K) .NE. 0) SH(K,I) = YMAX(I)
  205 CONTINUE
  206 CONTINUE
      RETURN
      END
      SUBROUTINE COTRA3
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      EQUIVALENCE (LIS1,  IFILE(7)),  (LIS2,   IFILE(8))
      EQUIVALENCE (ICRYS, IFILE(3)),  (ITS15, IFILE(15))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      PARAMETER (MXP=50, MXP1=MXP + 1)
      COMMON /PIEK/ RHO(MXP1), SH(3,MXP1), IBOTS(MXP1), BOTS(MXP1), NPK
      DIMENSION ISAME(MXP)
      COMMON /COTRA/ MULS(48,48), INDX(3,3), IDIMF, TT(3,48),
     *               EFST(48), FHST(48), PHST(48), CELLN(6), MAXHKL(3),
     *  P1SQ, PSQ, EMIN, EO2AV, SCX, BOVX, DAMPX, SMM, IOMAPS, SMAX,
     *  IRR(3,3,48), IORG, BBB, D2R, R2D,
     *  AMULT, ASYMM, ALATT, ASYMCL, NSYMC, ASYMC
      DIMENSION ORGVEC(3,8), SHD(3), SHDT(3), SHDTX(3), PLATT(3,4)
      LOGICAL SWORG
      DATA NORIG, SWORG /1 , .FALSE. /
      DO 238 I=1,NPK
      IF (IPOLA .EQ. 8) THEN
         SH(3,I) = - SH(2,I)
         SH(2,I) = SH(2,I) - SH(1,I)
         SH(3,I) = SH(3,I) - SH(1,I)
         SH(2,I) = SH(2,I) - SH(1,I)
         SH(1,I) = 0.0000001
         CALL KERNZA (0., PLATT, 3)
         CALL KERNZA (0., ORGVEC, 3)
         ENDIF
      DO 211 J=1,3
      SHJI = AMOD(2.0+SH(J,I),1.0)
  211 SH(J,I) = SHJI
      GOTO (212,220,215,230,230), IORG
  212 SH(1,I)=SH(1,I)/2.0
      SH(2,I)=SH(2,I)/2.0
  215 SH(3,I)=SH(3,I)/2.0
      GOTO 238
  220 IF (SH(1,I).LT.0.5) GOTO 215
      SH(1,I)=SH(1,I)-0.5
      SH2I = AMOD(SH(2,I)+0.5, 1.0)
      SH(2,I) = SH2I
      GOTO 215
  230 IF (SH(1,I).LT.0.5) GOTO 238
      DO 235 J=1,3
      SHJI = AMOD(SH(J,I)+0.5, 1.0)
      SH(J,I) = SHJI
  235 CONTINUE
  238 CONTINUE
      IF (IPOLA .EQ. 8) GOTO 242
      CALL RDCRYB (ICRYS, 'NORIG', KEND)
      IF (KEND .EQ. -1)
     *   CALL KERROR (' No NORIG record on CRYSDA file', 238, 'COTRA3')
      READ(CHIN, FMT = '(10X, I10)') NORIG
      DO 239 N = 1,NORIG
  239 READ (ICRYS, FMT = '(10X, 3F10.7)') (ORGVEC(J,N), J= 1,3)
      CALL FILCLO (ICRYS, 'KEEP')
      IF (NORIG * NLATT .EQ. 1) GOTO 248
      DO 241 N = 1, NLATT
      PLATT(1,N) = TLATT(1,N)
      PLATT(2,N) = TLATT(2,N)
      PLATT(3,N) = TLATT(3,N)
      IF (IPOLA.EQ.1 .OR. IPOLA.EQ.3 .OR. IPOLA.EQ.5) PLATT(1,N) = 0.
      IF (IPOLA.EQ.2 .OR. IPOLA.EQ.3 .OR. IPOLA.EQ.6) PLATT(2,N) = 0.
      IF (IPOLA.EQ.4 .OR. IPOLA.EQ.5 .OR. IPOLA.EQ.6) PLATT(3,N) = 0.
  241 CONTINUE
  242 DO 246 N1 = 1,NPK-1
      IF (RHO(N1) .LT. 0.0) GOTO 246
      DO 245 N2 = N1+1, NPK
      IF (RHO(N2) .LT. 0.0) GOTO 245
      DO 243 J = 1,3
  243 SHD(J) = SH(J,N1) - SH(J,N2)
      DO 244 IL = 1, NLATT
      SHDT(1) = SHD(1) + PLATT(1,IL)
      SHDT(2) = SHD(2) + PLATT(2,IL)
      SHDT(3) = SHD(3) + PLATT(3,IL)
      DO 244 N = 1,NORIG
      CALL DISTSQ (SHDT, ORGVEC(1,N), 0.4, SHDTX, DIST2)
      IF (DIST2 .LT. 0.16) THEN
         IF (.NOT. SWORG)
     *      WRITE (LIS2, FMT = '('' TEMP: Removed (close contact):'')')
         WRITE (LIS2,FMT='('' shiftvector'',I3,'' ='',3F8.4,'' ='',I3)')
     *      N2, (SH(J,N2),J=1,3), N1
         RHO(N2) = -9999.
         SWORG = .TRUE.
         GOTO 245
         ENDIF
  244 CONTINUE
  245 CONTINUE
  246 CONTINUE
      IF (.NOT. SWORG) GOTO 248
      NPK1 = 0
      DO 247 I = 1,NPK
      IF (RHO(I) .LT. 0.0) GOTO 247
      NPK1 = NPK1 + 1
      CALL KERNAB (SH(1,I), SH(1,NPK1), 3)
      RHO(NPK1) = RHO(I)
  247 CONTINUE
      NPK = NPK1
  248 CONTINUE
      CALL DISCHK (NA, ISAME, MXP)
      NBOTS = 0
      NBOTSP = 0
      DO 270 I=1,NPK
      IF (IBOTS(I) .NE. 0) NBOTS = NBOTS + 1
      IF (BOTS(I) .LE. 3.00) NBOTSP = NBOTSP + 1
  270 CONTINUE
      IF (NBOTS .EQ. 0) WRITE (LIS1, FMT='(/'' No close contacts''/)')
      IF (NBOTS .GT. 0) WRITE (LIS1, 275) NBOTS
  275 FORMAT (I4, ' of the shift vectors lead to bad contacts.'/
     *    '   Number of contacts (nr)  and shortests distance:'/
     *    42X, 'nr   dist')
      IF (NBOTS .EQ. 0 .AND. NBOTSP .GT. 0) WRITE (LIS1, FMT='(
     *    33X, ''Shortest distance:'')')
      RHOX = RHO(1) * .5
      WRITE(LIS1, FMT='('' Possible shift vectors:''/
     *  '' TR=   FOM    Tx       Ty       Tz'' )')
      NPKXX = NPK
      DMAX=2.999
      MSAME = NA / 3
      DO 280 I = 1, NPKXX
      IF (I .LE. 2 .OR. RHO(I) .GT. RHOX) NPK = I
      IF (IBOTS(I) .GT. 0) THEN
         IF (I .LE. 2 .OR. RHO(I) .GT. RHOX)
     *      WRITE (LIS1,276) I,RHO(I),(SH(J,I),J=1,3), IBOTS(I), BOTS(I)
  276    FORMAT (I4, F7.0, 3F9.5, I6, F7.2)
      ELSEIF (BOTS(I) .GT. DMAX) THEN
         IF (I .LE. 2 .OR. RHO(I) .GT. RHOX)
     *      WRITE (LIS1, 276) I, RHO(I), (SH(J,I),J=1,3)
      ELSE
         IF (I .LE. 2 .OR. RHO(I) .GT. RHOX)
     *      WRITE (LIS1, 277) I, RHO(I), (SH(J,I),J=1,3), BOTS(I)
  277    FORMAT (I4, F7.0, 3F9.5, 6X, F7.2)
         ENDIF
      IF (ISAME(I).GT.MSAME)
     *   WRITE (LIS1,FMT='('' Symmetry ?? shift nr. '', I3)') I
      IBOTS(I) = MIN0 ( IBOTS(I), 79)
  280 CONTINUE
      RHO(NPK+1)=-9999.0
      RETURN
      END
      SUBROUTINE R2TRAC (SEDEL2)
      DIMENSION SEDEL2(50,2)
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      EQUIVALENCE (IPR1, IFILE(6)), (LIS1, IFILE(7)), (LIS2, IFILE(8))
      EQUIVALENCE (IBINFO, IFILE(11)), (IBINFC, IFILE(12))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      PARAMETER (MAXBUF = 198)
      DIMENSION    BUFFC2(MAXBUF), BUFFO(MAXBUF),
     +             FITFC2(51),     FITFO(3)
      INTEGER      NITFC2, KENDFC, NITFO, KENDFO
      DIMENSION  BUF(2,24)
      EQUIVALENCE (SFTOT2,  FITFC2(2)),
     +            (SFPAR2, FITFC2(3)), (BUF(1,1),FITFC2(4))
      PARAMETER (MXP=50, MXP1=MXP + 1)
      COMMON /PIEK/ RHO(MXP1), X1(3,MXP1), IBOTS(MXP1), BOTS(MXP1), NPK
      COMMON /COTRA/ MULS(48,48), INDX(3,3), IDIMF, TT(3,48),
     *               EFST(48), FHST(48), PHST(48), CELLN(6), MAXHKL(3),
     *  P1SQ, PSQ, EMIN, EO2AV, SCX, BOVX, DAMPX, SMM, IOMAPS, SMAX,
     *  IRR(3,3,48), IORG, BBB, D2R, R2D,
     *  AMULT, ASYMM, ALATT, ASYMCL, NSYMC, ASYMC
      DIMENSION SHKL(3,24), SORTR2(50)
      WRITE (LIS1, 120)
  120 FORMAT (/' Calculate  R2  values for all  shift vectors'/
     *         ' (ignore possible overlap between molecules)')
      CALL BINIFF (1, IBINFO, 'BINFO', FITFO, NITFO, BUFFO, KENDFO)
      CALL BINIFF (1, IBINFC, 'BINFC2',FITFC2,NITFC2,BUFFC2,KENDFC)
      TOTREF = 0.
      SEOBS2 = 0.
      SEFS2 = 0.
      SEOBS4 = 0.
      CALL KERNZA (0.0, SEDEL2, 100)
      MPK = MIN0 (50, NPK)
  160 CALL BINIFF (0,IBINFC,'BINFC2',FITFC2,NITFC2,BUFFC2,KENDFC)
      IF (KENDFC .LT. 0) GOTO 500
      CALL BINIFF (0,IBINFO,'BINFO',FITFO,NITFO,BUFFO,KENDFO)
      IF (KENDFO .LT. 0) CALL KERNER(-2, 'R2TRAC')
      IF (FITFC2(1) .LT. 0.) GOTO 160
      CALL HKLC1U (FITFO(1), SHKL)
      CALL HKLSTL (SHKL, STL, STL2)
      IF (STL .GT. SMAX) GOTO 160
      FNORM = EXP(BOVX * STL2) / SQRT(SFTOT2 * ALATT)
      EF = SCX * FITFO(2) * FNORM
      IF (EF.GT.3.0) EF=3.0
      IF (EF.LT.EMIN) GOTO 160
      EF2 = EF**2
      SEOBS2 = SEOBS2 + EF2
      SEOBS4 = SEOBS4 + EF2**2
      CALL HKLEX1 (SHKL,  SHKL)
      DO 190 IS=1,NSYMM
      PHST(IS) = BUF(2,IS)
      IF (BUF(1,IS) .GE. -0.01) THEN
         EFST(IS) = BUF(1,IS)
      ELSE
         ISK = NINT (-BUF(1,IS))
         EFST(IS) = EFST(ISK)
         ENDIF
      PTRA = SHKL(1,1)*TT(1,IS)+ SHKL(2,1)*TT(2,IS)+ SHKL(3,1)*TT(3,IS)
      PHSTX    = PHST(IS) + 360. * PTRA
      PHST(IS) = AMOD (PHSTX, 360.)
  190 CONTINUE
      IF (ICENT.EQ.1) GOTO 196
      DO 195 IS = NSYMM+1 , NSYMC
      EFST(IS) = EFST(IS-NSYMM)
  195 PHST(IS) = 360. - PHST(IS-NSYMM)
  196 CONTINUE
      TOTREF = TOTREF + 1.
      DO 215 I = 1, MPK
      FA1 = 0.
      FB1 = 0.
      DO 210 IS=1,NSYMM
      PTRA =SHKL(1,IS)*X1(1,I)+ SHKL(2,IS)*X1(2,I)+ SHKL(3,IS)*X1(3,I)
      PTRA = 360.*PTRA + PHST(IS)
      PTRA = AMOD(PTRA, 360.)
      IF (PTRA .LT. 0.) PTRA = PTRA + 360.
                      FA1 = FA1 + EFST(IS) * COS(PTRA * D2R)
      IF (ICENT.EQ.1) FB1 = FB1 + EFST(IS) * SIN(PTRA * D2R)
  210 CONTINUE
      IF (ICENT.EQ.2) FA1 = FA1 + FA1
      EFS2 =(FA1**2 + FB1**2) * FNORM * FNORM
      IF (EFS2.GT.9.0) EFS2=9.0
      SEDEL2(I,1) = SEDEL2(I,1) + ( EF2 - EFS2 ) **2
      IF (I .EQ. 1) SEFS2 = SEFS2 + EFS2
  215 CONTINUE
      GOTO 160
  500 CONTINUE
      CALL FILCLO (IBINFC,'DELETE')
      CALL FILCLO (IBINFO,'KEEP')
      SEOBS2 = SEOBS2 / TOTREF
      SEFS2 = SEFS2 / TOTREF
      WRITE (LIS2, 910) SEOBS2,SEFS2
  910 FORMAT (/' Average of    observed' /
     +' EOBS**2     ',F10.4 /
     +' EFS**2      ',F10.4 ,' for first shift only')
      DO 919 I = 1, MPK
      SEDEL2(I,1) = SEDEL2(I,1) / SEOBS4
  919 CONTINUE
      WRITE (LIS1, 923) (I, SEDEL2(I,1), I=1, MPK)
      WRITE (LIS2, 923) (I, SEDEL2(I,1), I=1, MPK)
  923 FORMAT (/' R2 for all shift vectors '/
     *   '  nr  R2    nr  R2    ......'/
     *   5(I4, F6.3) )
      WRITE (LIS1, FMT='( '' R2 used for output selection.'' )')
      CALL KERNZA (0.0, SORTR2, 50)
      SORTR2(1) = SEDEL2(1,1)
      DO 933 I = 2, MPK
      DO 930 L = I, 2, -1
      IF (SEDEL2(I,1) .LT. SORTR2(L-1)) THEN
         SORTR2(L) = SORTR2(L-1)
         IF (L .EQ. 2) SORTR2(L-1) = SEDEL2(I,1)
      ELSE
         SORTR2(L) = SEDEL2(I,1)
         GOTO 933
         ENDIF
  930 CONTINUE
  933 CONTINUE
      WRITE (LIS2, 943) ( SORTR2(I), I=1, MPK)
  943 FORMAT (/' R2 sorted decreasing ... '/
     *   '      R2        R2        R2        R2        R2'/
     *   5(4X, F6.3) )
      R2MIN = SORTR2(1)
      XXX = ( 1. + R2MIN ) /2.
      IF (XXX .GT. 0.91) XXX = R2MIN * 1.10
      XXXR2 = XXX - R2MIN
      DO 950 I = 1, MPK
      SEDEL2(I,2) = ( XXX - SEDEL2(I,1) ) / XXXR2
      IF (SEDEL2(I,2) .LT. -9.99) SEDEL2(I,2) = -9.99
      IF (SEDEL2(I,2) .LE. 0.) GOTO 950
      IF (SEDEL2(I,1) .GT. 2. * R2MIN) SEDEL2(I,2) = -9.99
  950 CONTINUE
      IF (MPK .LE. 9) RETURN
      R2MAX = SORTR2(9)
      DO 953 I = 1, MPK
      IF (SEDEL2(I,1) .GT. R2MAX) SEDEL2(I,1) = 999. + SEDEL2(I,1)
  953 CONTINUE
      RETURN
      END
      SUBROUTINE SHIFTR (SEDEL2)
      DIMENSION SEDEL2(50,2)
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      LOGICAL SWPRI
      EQUIVALENCE (SWPRI, SWITCH(10))
      EQUIVALENCE (IATOMS, IFILE(2)), (IATOLD, IFILE(10))
      EQUIVALENCE (IPR1, IFILE(6)), (LIS1, IFILE(7)), (LIS2, IFILE(8))
      PARAMETER (KUSER2=30000)
      PARAMETER (MAXAT=993)
      COMMON / / XX(KUSER2), ATXYZ(10,MAXAT), NATX
      DIMENSION BLACOM(42000)
      EQUIVALENCE (BLACOM(1), XX(1))
      COMMON /ATNAMA/ ATNAME(MAXAT)
      CHARACTER * 6  ATNAME
      PARAMETER (MXP=50, MXP1=MXP + 1)
      COMMON /PIEK/ RHO(MXP1), SH(3,MXP1), IBOTS(MXP1), BOTS(MXP1), NPK
      DIMENSION XLOCK(3)
      DATA DMAXTR / 0.431 /
      IF (SWPRI) WRITE (LIS2, FMT = ' ('' The new coordinates are: '' /
     *      '' Atom name      x         y         z'' /)')
      DO 123 N = 1,NATX
      ATXYZ(8,N) = ATXYZ(1,N)
      ATXYZ(9,N) = ATXYZ(2,N)
  123 ATXYZ(10,N) = ATXYZ(3,N)
      REWIND IATOMS
      NOTR2 = 0
      DO 540 I = 1, NPK
      INEW = I - NOTR2
      DO 330 N = 1,NATX
      ATXYZ(1,N) = ATXYZ(8,N) + SH(1,I)
      ATXYZ(2,N) = ATXYZ(9,N) + SH(2,I)
      ATXYZ(3,N) = ATXYZ(10,N) + SH(3,I)
      CALL LOCKIN (ATXYZ(1,N), DMAXTR, XLOCK, DISTTR, NPOSTR)
      IF (NPOSTR .GT. 1) THEN
         IF (I .EQ. 1) WRITE (LIS1, 226) N, ATNAME(N)
  226    FORMAT (' Shift nr 1, Atom', I4,' = ', A6,
     *      ' lies on a special position')
         CALL KERNAB (XLOCK, ATXYZ(1,N), 3)
         IF (I .EQ. 1 .AND. DISTTR .GT. 0.2) WRITE (LIS1, FMT ='
     *      ('' Atom locked-in: dist ='', F5.2, '' ang.'')') DISTTR
         IF (DISTTR .GT. 0.2) WRITE (LIS2, FMT ='
     *      ('' Atom locked-in: dist ='', F5.2, '' ang.'')') DISTTR
         ENDIF
      IF (SWPRI) WRITE (LIS2, FMT = ' (4X, A6, 3F10.5)' )
     *      ATNAME(N), (ATXYZ(J,N), J=1,3)
  330 CONTINUE
      R299 = SEDEL2(I,1)
      IF (R299 .GT. 999.) R299 = R299 - 999.
      IF (SEDEL2(I,2) .LE. 0.) THEN
         NOTR2 = NOTR2 + 1
         WRITE (LIS1, 332) I, R299, SEDEL2(I,2)
         GOTO 540
      ELSEIF (SEDEL2(I,1) .LT. 99.) THEN
         WRITE (LIS1, 331) I, R299, SEDEL2(I,2), INEW
  331    FORMAT ( ' TRACOR shift #',I2, '  R2 =', F6.3,
     *            '   Q2=',F5.2, '   accepted:  TR=', I2)
      ELSEIF (I .LT. 9) THEN
         WRITE (LIS1, 1332) I, R299, SEDEL2(I,2), INEW
 1332    FORMAT ( ' TRACOR shift #',I2, '  R2 =', F6.3,
     *       '   Q2=',F5.2, '   accepted:  TR=',I2, '  R2=BIG !')
      ELSE
         NOTR2 = NOTR2 + 1
         WRITE (LIS1, 332) I, R299, SEDEL2(I,2)
  332    FORMAT ( ' TRACOR shift #',I2, '  R2 =', F6.3,
     *       '   Q2=',F5.2, '   rejected         R2 is too BIG')
         GOTO 540
         ENDIF
      WRITE (CHOUT, 333) INEW, (SH(J,I),J=1,3),
     *   SEDEL2(I,1), RHO(I), IBOTS(I)
  333 FORMAT ('TR=', I3, ' Shift:', 3F8.4, ' R2=', F6.3,
     *       ' FOM=', F6.0, ' X= ', I2)
      IFOM = NINT(RHO(I))
      IF (SEDEL2(I,1) .GT. 99.) SEDEL2(I,1) = SEDEL2(I,1) - 999.
      CALL ATOMWT (IATOMS, ATXYZ, ATNAME, NATX, INEW,
     *   SEDEL2(I,1), IFOM, IBOTS(I))
  540 CONTINUE
      INEW = NPK - NOTR2
      CALL FILCLO (IATOMS, 'KEEP')
      WRITE (IPR1, 544) NPK
  544 FORMAT (' Number of shift vectors found by TRACOR :', I3)
      IF (INEW .LT. NPK) WRITE (IPR1, 545) INEW
  545 FORMAT (' Number of vectors accepted after R2 test:', I3)
      IF (IBOTS(1) .GT. 0) THEN
         WRITE (IPR1, 550)
         WRITE (LIS1, 550)
  550    FORMAT(/' NOTE: the first shift vector brings the fragment in'
     *         /'       collision with a symmetry related molecule!'/)
         MSAME = NATX / 3
         IF (IBOTS(1) .GT. MSAME) THEN
            WRITE (IPR1, 554)
            WRITE (LIS1, 554)
  554      FORMAT(' The fragment is probably put on a symmetry element:'
     *           /'       check if its internal symmetry is correct.'/)
            ENDIF
         ENDIF
      IF (INEW .GT. 1) WRITE (LIS1, 120)
  120 FORMAT (/' All shift vectors have been applied and tested. '/
     * ' All accepted parameter sets are written to the ATOMS file,' /
     * ' and transferred to TRAVEC'/
     * ' but TRAVEC is going to reorder the atom sets.'/)
      RETURN
      END
      SUBROUTINE ATOMWT (IATOMS, ATXYZ,ATNAME,NAT, NNNN, R2,IFOM,IBOTS)
      DIMENSION ATXYZ(10,NAT), ATNAME(NAT)
      CHARACTER *6 ATNAME
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      EQUIVALENCE    (LIS2, IFILE(8))
      CHARACTER *6 CBOTS, TBOTS
      IF (NNNN .EQ. 1) WRITE (LIS2, FMT='(/
     *   '' Output ATOMS    header lines''/)')
      IF (IBOTS .GT. 0) THEN
         CALL KERI2C (IBOTS, TBOTS, 2)
         CBOTS(5:6) = TBOTS(1:2)
         CBOTS(1:4) = ' X= '
      ELSE
         CBOTS = ' '
         ENDIF
      IT = MAX0 (0, 10000 * (ITIME(1)-1900) + 100 * ITIME(2) + ITIME(3))
      WRITE (IATOMS, 102) CCODE, IT, KEYS(13), NNNN, R2, IFOM, CBOTS
  102 FORMAT ('ATOMS ', A6, ' < TRACOR :'
     *  , I7, ' RUN', I4, ' TR=', I3, ' R2= ', F5.3, ' FOM=', I5, A6)
      WRITE (LIS2, 1102) CCODE, IT, KEYS(13), NNNN, R2, IFOM, CBOTS
 1102 FORMAT (/' ', 'ATOMS ', A6, ' < TRACOR :'
     *  , I7, ' RUN', I4, ' TR=', I3, ' R2= ', F5.3, ' FOM=', I5, A6)
      IF (CHOUT .NE. ' ') THEN
         WRITE (IATOMS, FMT = '(''REMARK '', A65)') CHOUT(1:65)
         WRITE (LIS2, FMT = '('' REMARK '', A65)') CHOUT(1:65)
         CHOUT = ' '
         ENDIF
      DO 109 I = 1, NAT
      WRITE (IATOMS, 104) ATNAME(I), (ATXYZ(J,I), J=1,3)
  104 FORMAT ('ATOM  ', A6, 2X, 3F9.5)
  109 CONTINUE
      WRITE (IATOMS, FMT = '(''END'')')
      RETURN
      END
      SUBROUTINE D1FOUR (IO,NGR1,SCALE)
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      EQUIVALENCE (LIS2,  IFILE( 8))
      EQUIVALENCE (ITS15, IFILE(15))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      PARAMETER (MNGR1=200, MNGR2=163, NPKM1=50, NPKM2=100, LC24=24)
      COMMON / / SI(1250),  RHOSUM(MNGR1),   JUNK(MNGR2*MNGR2-MNGR1),
     +           RHOMAX(NPKM2), IC1(MNGR1), IC2(MNGR1),  IMAX(NPKM2),
     +           JMAX(NPKM2), XMAX(NPKM2), YMAX(NPKM2),  RHOP(NPKM2),
     +           IRHO(LC24)
      DIMENSION BLACOM(42000)
      EQUIVALENCE (BLACOM(1), SI(1))
      PARAMETER (MXP=50, MXP1=MXP + 1)
      COMMON /PIEK/ RHO(MXP1), TX(3,MXP1), IBOTS(MXP1), BOTS(MXP1), NPK
      DIMENSION CO(1000)
      EQUIVALENCE (SI(251),CO(1))
      CALL KERNZA (0.,RHOSUM,MNGR1)
      PA = 1000. / FLOAT(NGR1)
      DO 40 I1=1,NGR1
   40 IC1(I1) = FLOAT(I1) * PA + 0.5
   60 READ(ITS15) IH,FA,FB
      IF (IH.EQ.999) GOTO 70
      DO 50 I1=1,NGR1
      II = MOD (IH * IC1(I1), 1000)
      IF (II .LE. 0) II = II + 1000
      RHOSUM(I1)=RHOSUM(I1)+FA*CO(II)+FB*SI(II)
   50 CONTINUE
      GOTO 60
   70 CALL KERNZA (0.,RHOMAX,NPKM1)
      CALL KERNZI (0,IMAX,NPKM1)
      SCALE=0.0
      DO 100 I1=1,NGR1
      SCALE=SCALE+RHOSUM(I1)**2
      DO 110 I3=1,NPKM1
      IF (RHOSUM(I1).GT.RHOMAX(I3)) GOTO 130
  110 CONTINUE
      GOTO 100
  130 DO 120 I4=1,NPKM1
      IF ((NPKM1+1-I4).EQ.I3) GOTO 210
      IMAX(NPKM1+1-I4)=IMAX(NPKM1-I4)
      RHOMAX(NPKM1+1-I4)=RHOMAX(NPKM1-I4)
  120 CONTINUE
  210 IMAX(I3)=I1
      RHOMAX(I3)=RHOSUM(I1)
  100 CONTINUE
      DO 170 I1=1,NPKM1-1
      IF(RHOMAX(I1).LE. 0.01) GOTO 170
      I4=I1+1
      DO 180 I2=I4,NPKM1
      IF (RHOMAX(I2).LE. .01) GOTO 180
      IDIF1=IABS(IMAX(I1)-IMAX(I2))
      IDIF1=MIN0(IDIF1,NGR1-IDIF1)
      IF (IDIF1.LT.3 ) GOTO 200
      GOTO 180
  200 RHOMAX(I2)=0.
  180 CONTINUE
  170 CONTINUE
      SCALE = 100. * SQRT (FLOAT(NGR1) / SCALE)
      IF (IO.LE.0) GOTO 330
      WRITE (LIS2,1003)
 1003 FORMAT ('      ***1D.-FOURIER MAP***  '/)
      ITW2=1
      ITW =1
      WRITE (LIS2,1008)
 1008 FORMAT('  U= 1 ,   2 ,   3....ETC'/)
      GOTO 310
  340 ITW2=ITW2+1
      ITW=LC24*(ITW2-1)+1
  310 NPAG=LC24*ITW2
      IF (NGR1.LT.NPAG) NPAG=NGR1
      I3=0
      DO 1010 I1=ITW,NPAG
      I3=I3+1
      IRHO(I3)=IFIX(RHOSUM(I1)*SCALE+.5)
 1010 CONTINUE
      WRITE (LIS2,1000)  (IRHO(I1),I1=1,I3)
 1000 FORMAT (1X,24I5)
      IF (NGR1.GT.NPAG) GOTO 340
  330 N=0
      DO 250 I1=1,NPKM1
      IF (RHOMAX(I1).LE. .01) GOTO 250
      IX1=IMAX(I1)-1
      IX2=IMAX(I1)
      IX3=IMAX(I1)+1
      IXX3=IX3
      IXX1=IX1
      IF (IX2.EQ.NGR1) IXX3=IX3-NGR1
      IF (IX2.EQ.1) IXX1=IX1+NGR1
      RHO1=RHOSUM(IXX1)
      RHO2=RHOMAX(I1)
      RHO3=RHOSUM(IXX3)
      IF (RHO1.GT.RHO2 .OR. RHO3.GT.RHO2) GOTO 250
      CALL INTPOL(RHO1,RHO2,RHO3,IX1,IX3,X,NGR1)
      XMAX(I1)=AMOD(X,1.0)
      RHOP(I1)=RHO2*SCALE
      N=N+1
      RHOP(N)=RHOP(I1)
      XMAX(N)=XMAX(I1)
      IF (N .EQ. MXP) GOTO 260
  250 CONTINUE
  260 RHOP(N+1)=-9999.
      XMAX(N+1) = 0.
      NPK = N
      RETURN
      END
      SUBROUTINE D2FOUR (IO, NGR1, NGR2, SCALE)
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      EQUIVALENCE (LIS2,  IFILE( 8))
      EQUIVALENCE (ITS15, IFILE(15))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      PARAMETER (MNGR1=200, MNGR2=163, NPKM2=100, LC24=24)
      COMMON / / SI(1250),RHOSUM(MNGR2,MNGR2),RHOMAX(NPKM2),IC1(MNGR1),
     +           IC2(MNGR1),  IMAX(NPKM2),  JMAX(NPKM2),   XMAX(NPKM2),
     +           YMAX(NPKM2),RHOP(NPKM2),IRHO(LC24)
      DIMENSION BLACOM(42000)
      EQUIVALENCE (BLACOM(1), SI(1))
      PARAMETER (MXP=50, MXP1=MXP + 1)
      COMMON /PIEK/ RHO(MXP1), TX(3,MXP1), IBOTS(MXP1), BOTS(MXP1), NPK
      DIMENSION CO(1000),AH(MNGR1),BH(MNGR1)
      EQUIVALENCE (SI(251),CO(1)),(AH(1),IMAX(1)),(BH(1),XMAX(1))
      CALL KERNZA (0.,RHOSUM,MNGR2*MNGR2)
      PA1 = 1000. / FLOAT(NGR1)
      PA2 = 1000. / FLOAT(NGR2)
      DO 40 I1=1,NGR1
   40 IC1(I1) = FLOAT(I1) * PA1 + 0.5
      DO 42 I2=1,NGR2
   42 IC2(I2) = FLOAT(I2) * PA2 + 0.5
      IHO=-1
   48 READ(ITS15) IH,IK,FA,FB
      IF (IHO.EQ.-1) GOTO 53
      IF (IH.EQ.IHO) GOTO 54
      DO 52 I1=1,NGR1
      IHX = MOD (IHO * IC1(I1), 1000)
      IF (IHX .LE. 0) IHX = IHX + 1000
      DO 52 I2=1,NGR2
   52 RHOSUM(I1,I2)=RHOSUM(I1,I2)+AH(I2)*CO(IHX)+BH(I2)*SI(IHX)
      IF (IH.EQ.999) GOTO 70
   53 CALL KERNZA(0.,AH,MNGR1)
      CALL KERNZA(0.,BH,MNGR1)
      IHO=IH
   54 DO 55 I2=1,NGR2
      KY = MOD (IK * IC2(I2), 1000)
      IF (KY .LE. 0) KY = KY + 1000
      AH(I2)=AH(I2) + FA*CO(KY) + FB*SI(KY)
   55 BH(I2)=BH(I2) + FB*CO(KY) - FA*SI(KY)
      GOTO 48
   70 CONTINUE
      CALL KERNZA (0.,RHOMAX,NPKM2)
      CALL KERNZI (0,IMAX,NPKM2)
      CALL KERNZI (0,JMAX,NPKM2)
      SCALE=0.0
      DO 100 I1=1,NGR1
      DO 100 I2=1,NGR2
      SCALE=SCALE+RHOSUM(I1,I2)**2
      IF (RHOSUM(I1,I2).LE.RHOMAX(NPKM2)) GOTO 100
      DO 110 I3=1,NPKM2
      IF (RHOSUM(I1,I2).GT.RHOMAX(I3)) GOTO 130
  110 CONTINUE
      GOTO 100
  130 DO 120 I4=1,NPKM2
      IF ((NPKM2+1-I4).EQ.I3) GOTO 210
      IMAX(NPKM2+1-I4)=IMAX(NPKM2-I4)
      JMAX(NPKM2+1-I4)=JMAX(NPKM2-I4)
      RHOMAX(NPKM2+1-I4)=RHOMAX(NPKM2-I4)
  120 CONTINUE
  210 IMAX(I3)=I1
      JMAX(I3)=I2
      RHOMAX(I3)=RHOSUM(I1,I2)
  100 CONTINUE
      DO 170 I1=1,NPKM2-1
      IF (RHOMAX(I1) .LE. .01) GOTO 170
      I4=I1+1
      DO 180 I2=I4,NPKM2
      IF (RHOMAX(I2) .LE. .01) GOTO 180
      IDIF1=IABS(IMAX(I1)-IMAX(I2))
      IDIF1=MIN0(IDIF1,NGR1-IDIF1)
      IDIF2=IABS(JMAX(I1)-JMAX(I2))
      IDIF2=MIN0(IDIF2,NGR2-IDIF2)
      IF (IDIF1.GE.3 .OR.  IDIF2.GE.3) GOTO 180
      RHOMAX(I2)=0.
  180 CONTINUE
  170 CONTINUE
      SCALE = 100. * SQRT (FLOAT(NGR1 * NGR2) / SCALE)
      IF (IO.LE.0) GOTO 330
      ITW2=1
      ITW=1
      GOTO 310
  340 ITW2=ITW2+1
      ITW=LC24*(ITW2-1)+1
  310 NPAG=LC24*ITW2
      IF (NGR1.LT.NPAG) NPAG=NGR1
      WRITE (LIS2,1012) (I,I=ITW,NPAG)
 1012 FORMAT ('1',3X,'U=',24(2X,I3))
      WRITE (LIS2,1013)
 1013 FORMAT(' ','V=')
      DO 300 I1=1,NGR2
      I3=0
      DO 302 I2=ITW,NPAG
      I3=I3+1
      IRHO(I3)=IFIX(RHOSUM(I2,I1)*SCALE+.5)
  302 CONTINUE
      WRITE (LIS2,1000) I1,(IRHO(I2),I2=1,I3)
 1000 FORMAT ('0',I3,2X,24I5)
  300 CONTINUE
      IF (NGR1.GT.NPAG) GOTO 340
  330 N=0
      DO 250 I1=1,NPKM2
      IF (RHOMAX(I1).LE. .01) GOTO 250
      IX1=IMAX(I1)-1
      IX2=IMAX(I1)
      IX3=IMAX(I1)+1
      IY1=JMAX(I1)-1
      IY2=JMAX(I1)
      IY3=JMAX(I1)+1
      IXX3=IX3
      IXX1=IX1
      IF (IX2.EQ.NGR1) IXX3=IX3-NGR1
      IF (IX2.EQ.1) IXX1=IX1+NGR1
      RHO1=RHOSUM(IXX1,IY2)
      RHO2=RHOMAX(I1)
      RHO3=RHOSUM(IXX3,IY2)
      IF (RHO1.GT.RHO2 .OR. RHO3.GT.RHO2) GOTO 250
      CALL INTPOL(RHO1,RHO2,RHO3,IX1,IX3,X,NGR1)
      IM1=IMAX(I1)
      XMAX(I1)=AMOD(X,1.0)
      RHOP(I1)=RHO2
      IYY1=IY1
      IYY3=IY3
      IF (IY2.EQ.NGR2) IYY3=IY3-NGR2
      IF (IY2.EQ.1)    IYY1=IY1+NGR2
      RHO1=RHOSUM(IX2,IYY1)
      RHO2=RHOMAX(I1)
      RHO3=RHOSUM(IX2,IYY3)
      IF (RHO1.GT.RHO2 .OR. RHO3.GT.RHO2) GOTO 250
      CALL INTPOL (RHO1,RHO2,RHO3,IY1,IY3,X,NGR2)
      YMAX(I1)=AMOD(X,1.0)
      RHOP(I1)=(RHOP(I1)+RHO2)/2.*SCALE
      N=N+1
      RHOP(N)=RHOP(I1)
      XMAX(N)=XMAX(I1)
      YMAX(N)=YMAX(I1)
      IF (N .EQ. MXP) GOTO 260
  250 CONTINUE
  260 RHOP(N+1)=-999.
      XMAX(N+1) = 0.
      YMAX(N+1) = 0.
      NPK = N
      IF (NPK .LE. 1) RETURN
      IF (ILAUE .NE. 8) RETURN
      RHM = 0.45 * RHOP(1)
      NPK = 0
      DO 550 I1 = 1, N
      IF (RHOP(I1) .LT. RHM) GOTO 550
      NPK = NPK + 1
      RHOP(NPK) = RHOP(I1)
      XMAX(NPK) = XMAX(I1)
      YMAX(NPK) = YMAX(I1)
      IF (I1 .EQ. N) GOTO 550
      RHM1 = 0.83 * RHOP(I1)
      DO 540 I2 = I1+1, N
      IF (RHOP(I2) .LT. RHM1) GOTO 540
      XM1 = XMAX(I1) - XMAX(I2)
      IF (XM1 .LT. 0.0) XM1 = XM1 + 1.0
      YM1 = YMAX(I1) - YMAX(I2)
      IF (YM1 .LT. 0.0) YM1 = YM1 + 1.0
      IF (ABS(XM1-0.333).LE.0.02 .AND. ABS(YM1-0.667).LE.0.02) GOTO 530
      IF (ABS(XM1-0.667).GT.0.02 .OR.  ABS(YM1-0.333).GT.0.02) GOTO 540
  530 RHOP(I2) = 0.0
  540 CONTINUE
  550 CONTINUE
      RETURN
      END
      SUBROUTINE INTPOL(RHO1,RHO2,RHO3,IX1,IX3,X,IGRID)
      S=1.
      IF (RHO1.LT.RHO3) GOTO 10
      RHO33=RHO1-RHO3
      RHO2=RHO2-RHO3
      RHO11=0.
      IX1=IX3
      S=-1.
      GOTO 20
   10 RHO33=RHO3-RHO1
      RHO2=RHO2-RHO1
      RHO11=0.
   20 TEL=(RHO33/RHO2)-4.
      A=(.5*RHO33)-RHO2
      B=RHO2-A
      DENOM=(2.*RHO33/RHO2)-4.
      X=TEL/DENOM
      RHOM = (A*X+B)*X
      IF (NINT(S) .EQ. -1) GOTO 40
      RHOM=RHOM+RHO1
      GOTO 50
   40 RHOM=RHOM+RHO3
   50 X=(FLOAT(IX1)+S*X)/FLOAT(IGRID)
      RHO2=RHOM
      RETURN
      END
      SUBROUTINE D3FOUR
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      EQUIVALENCE (ICRYS, IFILE(3)), (ICON, IFILE(4))
      EQUIVALENCE (IBINBI, IFILE(18))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      PARAMETER (KUSER2=30000, KUSER1=KUSER2/3)
      PARAMETER (MXP=50, MXP1=MXP + 1)
      COMMON /PIEK/ RHO(MXP1),  X(3,MXP1), IBOTS(MXP1), BOTS(MXP1), NPK
      COMMON /COTRA/ MULS(48,48), INDX(3,3), IDIMF, TT(3,48),
     *               EFST(48), FHST(48), PHST(48), CELLN(6), MAXHKL(3),
     *  P1SQ, PSQ, EMIN, EO2AV, SCX, BOVX, DAMPX, SMM, IOMAPS, SMAX,
     *  IRR(3,3,48), IORG, BBB, D2R, R2D,
     *  AMULT, ASYMM, ALATT, ASYMCL, NSYMC, ASYMC
      COMMON /FFTDA/ SCALE, MH(3), NPP(3), XLMIN(3), XLMAX(3)
      CALL KERNAI (MAXHKL, MH, 3)
      CALL KERNAB (CELLN, CELL, 6)
      NSYMM = 1
      IF (ILAUE .EQ. 5) THEN
      NLATT = 2
      IMULT = 2
      TLATT(1,2) = 0.5
      TLATT(2,2) = 0.5
      TLATT(3,2) = 0.
      ELSE
      NLATT = 1
      IMULT = 1
      ENDIF
      ICENT = 1
      ILAUE = 1
      ISYST = 1
      ILATT = 1
      CALL RCELLR (CELL, VOLUM, RCELL)
      CALL CELLRR (CELL, RRMAT)
      CALL MATF2C (CELL, FRAC2C)
      KEYS(27) = 5
      CALL FFTRIN (KUSER1)
      CALL PP1
      CALL SEART
      CALL RDCRYS (ICRYS)
      RETURN
      END
      SUBROUTINE FFTRIN (KUSER1)
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      EQUIVALENCE (ICRYS, IFILE(3)), (ICON, IFILE(4))
      EQUIVALENCE (LIS1, IFILE(7)), (LIS2, IFILE(8))
      EQUIVALENCE (IBINFF, IFILE(16)), (IFMAP,  IFILE(17))
      EQUIVALENCE (ISCRA,  IFILE(18))
      EQUIVALENCE (KEYS(27), IMAP), (KEYS(28), IHALF)
      LOGICAL      SWPRI, PRIMAP
      EQUIVALENCE (SWPRI, SWITCH(10)), (PRIMAP, SWITCH(11))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     *               WAVE,     CELALL(10),  AMOLW,      ZET,
     *               NELEC,    F000,        ABSMU,      ICENT,
     *               ILATT,    ISYST,       ILAUE,      IMULT,
     *               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     *         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     *         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      COMMON /CRYSB/ SPGR,     WAVEAT,      CELATY(10)
      CHARACTER      SPGR *16, WAVEAT *2,   CELATY *2
      COMMON /FFTDA/ SCALE, MH(3), NPP(3), XLMIN(3), XLMAX(3)
      COMMON /SEARTA/ D2R, DMPIC, DMAXB, DMOUT, DMINB, ANGM(2), MCON,
     *        SEARDX, NPIC, NATIN, NAT, NATS, NRECY, PSQ, NATREC, KPROG
      PARAMETER (MAXBUF = 198)
      DIMENSION FITFFT(5), BUFFFT(MAXBUF), IGM3(3)
      DIMENSION MAXHKL(3)
      FACTOR  = 0.3
      CALL KERNZA (0.0, XLMIN, 3)
      CALL KERNZA (1.0, XLMAX, 3)
      IGM3(1) = 1
      IGM3(2) = 1
      IGM3(3) = 2
      CALL BINIFF (1, IBINFF, 'BINFFT', FITFFT, NITFFT, BUFFFT, NEND)
      PSQ = BUFFFT(27)
      CALL FILINQ (ISCRA, 'BINBIG', 'UNFORMATTED', 'OUTPUT', KINQ)
      CALL FILINQ (IFMAP, 'FMAPT', 'UNFORMATTED', 'OUTPUT', KINQ)
      CALL KERF2I (BUFFFT(22), MAXHKL, 3)
      IF (MAXHKL(1) .EQ. 0) CALL KERF2I (BUFFFT(7), MAXHKL, 3)
      CALL KERNAI (MH, MAXHKL, 3)
      DO 105 I=1,3
      MH(I) = MIN0 (MH(I), MAXHKL(I))
  105 IF (MH(I) .LE. 0) MH(I) = 1
      IHALF = 0
      IHALF = 1
      IF (ILATT .EQ. 1 .OR. ILATT .EQ. 3)  IHALF = 0
      FACTOR = 999.
      DO 219 I=1,3
      NPP(I) = 2 * MH(I) + 2
  219 FACTOR = AMIN1 (FACTOR, 0.90 * CELL(I)/NPP(I))
      IF (IHALF.NE.0 .AND. IGM3(3).EQ.1) IGM3(2) = 2
      GOTO 242
  240 DO 241 I=1,3
  241 NPP(I) = CELL(I) / FACTOR + 0.5
  242 DO 280 I=1,3
      ISGG = MOD (NPP(I), IGM3(I))
      IF (ISGG.NE.0) NPP(I) = NPP(I) + IGM3(I) - ISGG
  250 NTEST = NPP(I)
      DO 270 J=2,5
  260 IF (NTEST.NE.(NTEST/J)*J) GOTO 270
      NTEST = NTEST / J
      IF (NTEST.EQ.1) GOTO 280
      GOTO 260
  270 CONTINUE
      NPP(I) = NPP(I) + IGM3(I)
      GOTO 250
  280 CONTINUE
      WRITE (LIS2, 282) FACTOR
  282 FORMAT (' The GRID spacing is approximately', F6.3, ' Angstrom')
      IF (NPP(1) .LE. 250) GOTO 400
      WRITE (LIS2, 320)
  320 FORMAT (/' NX GREATER THAN 250 (SEE SUBR. -OUTPUT-). RESET.'/)
      FACTOR = FACTOR * FLOAT(NPP(1)) / 245.
      GOTO 240
  400 I = (NPP(1)+2) * (NPP(3)+2)
      IF (I.LT.KUSER1) GOTO 406
      FACTOR = FACTOR * 1.02 * SQRT(FLOAT(I)/FLOAT(KUSER1))
      WRITE (LIS1, 405)
      WRITE (LIS2, 405)
  405 FORMAT (' TOO MANY GRID POINTS FOR PEAK SEARCH. RESET.'/)
      GOTO 240
  406 WRITE (LIS2, 407) MH, NPP
  407 FORMAT (/
     + ' Maximum indices allowed   h:', I4, '    k:', I4, '    l:', I4/
     + ' Number of grid points    Nx:', I4, '   Ny:', I4, '   Nz:', I4)
      DO 408 I = 1,3
      IF (NPP(I) .GE. 2 * MH(I) + 2) GOTO 408
      MH(I) = NPP(I) / 2 -1
      WRITE (LIS2, FMT='('' Reset MAXHKL:'')')
      GOTO 406
  408 CONTINUE
      IF (SWPRI .AND. PRIMAP) WRITE (LIS2, 410) XLMAX
  410 FORMAT (' FOURIER MAP TO BE PRINTED FROM:'/' X =  0.0  TO',
     +         F7.3,',  Y =  0.0  TO',F7.3,',  Z =  0.0  TO',F7.3)
      RETURN
      END
          SUBROUTINE SEART
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      EQUIVALENCE (ICON,  IFILE( 4))
      EQUIVALENCE (LIS1,  IFILE( 7)), (LIS2, IFILE(8))
      EQUIVALENCE (IFMAP, IFILE(17))
      LOGICAL SWPRI
      EQUIVALENCE (SWPRI, SWITCH(10))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      PARAMETER (KUSER2=30000)
      PARAMETER (MAXAT=993)
      COMMON / / NSTORE(KUSER2), X(4,MAXAT)
      DIMENSION NBLACO(42000)
      EQUIVALENCE (NBLACO(1), NSTORE(1))
      PARAMETER (MXP=50, MXP1=MXP + 1)
      COMMON /PIEK/ RHO(MXP1), XT(3,MXP1), IBOTS(MXP1), BOTS(MXP1), NPK
      COMMON /SEARTA/ D2R, DMPIC, DMAXB, DMOUT, DMINB, ANGM(2), MCON,
     *        SEARDX, NPIC, NATIN, NAT, NATS, NRECY, PSQ, NATREC, KPROG
      DIMENSION ITLE(20)
      D2R = ATAN(1.0) / 45.0
      DMPIC = 0.5
      REWIND IFMAP
      READ (IFMAP) ITLE,IMAP,IHALF
      WRITE (LIS1, 192)
  192 FORMAT(/' FFT PEAK SEARCH FOR TRACOR === VERSION JUN 1987')
      NPIC = MXP
      NAT  = MXP
      CALL PKSRT (X, MAXAT, IHALF, IFMAP)
      CALL FILCLO (IFMAP, 'DELETE')
      IF (SWPRI) WRITE (LIS2, 445)
  445 FORMAT ('0Untransformed peaks from 3D Fourier map'/
     *     '   PEAK',7X,'HEIGHT',13X,'X',9X,'Y',9X,'Z')
      RHOMIN = 0.50 * X(4,1)
      DO 520 I=1,NPIC
      IF (X(4,I) .LT. RHOMIN) GOTO 530
      DO 515 J=1,3
  515 XT(J,I)=X(J,I)
      RHO(I)=X(4,I)
      IF (SWPRI) WRITE (LIS2, 517) I,RHO(I),(XT(J,I),J=1,3)
  517 FORMAT (' ' ,3X,I2,7X,F6.0,10X,F7.4,3X,F7.4,3X,F7.4)
      IF (I .EQ. 50) GOTO 550
  520 CONTINUE
      I = NPIC + 1
  530 I = I - 1
  550 DO 555 J=1,3
  555 XT(J,I+1)=0.
      RHO(I+1)=-9999.
      NPK = I
      RETURN
      END
      SUBROUTINE DISCHK (NA, ISAME, MXP2)
      DIMENSION ISAME(MXP2)
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      EQUIVALENCE (IATOMS, IFILE(2))
      EQUIVALENCE (LIS1, IFILE(7)), (LIS2, IFILE(8))
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      PARAMETER (KUSER2=30000, KUSER1=KUSER2/3)
      PARAMETER (MAXAT=993)
      COMMON / / XX(KUSER2), ATXYZ(10,MAXAT), NATX
      DIMENSION IZAT(MAXAT), WORK(3,KUSER1)
      EQUIVALENCE (WORK(1,1), IZAT(1), XX(1))
      DIMENSION BLACOM(42000)
      EQUIVALENCE (BLACOM(1), XX(1))
      COMMON /ATNAMA/ ATNAME(MAXAT)
      CHARACTER * 6  ATNAME
      PARAMETER (MXP=50, MXP1=MXP + 1)
      COMMON /PIEK/ RHO(MXP1), SH(3,MXP1), IBOTS(MXP1), BOTS(MXP1), NPK
      COMMON /COTRA/ MULS(48,48), INDX(3,3), IDIMF, TT(3,48),
     *               EFST(48), FHST(48), PHST(48), CELLN(6), MAXHKL(3),
     *  P1SQ, PSQ, EMIN, EO2AV, SCX, BOVX, DAMPX, SMM, IOMAPS, SMAX,
     *  IRR(3,3,48), IORG, BBB, D2R, R2D,
     *  AMULT, ASYMM, ALATT, ASYMCL, NSYMC, ASYMC
      DIMENSION TL(3,4)
      DIMENSION COSA(3), SINA(3), Q(3,3), SSH(3), XY(3), XGRAV(3)
      DIMENSION RATIO(MXP)
      CALL KERNZA(0.0, TL, 12)
      GOTO (37,31,31,31,31,34,33),ILATT
   31 DO 32 I=1,3
   32 TL(I,2)=0.5
      IF(ILATT.LT.5) TL(ILATT-1,2)=0.0
      GOTO 37
   33 TL(1,2)=1.0/3.0
      TL(2,2)=2*TL(1,2)
      TL(3,2)=TL(2,2)
      TL(1,3)=TL(2,2)
      TL(2,3)=TL(1,2)
      TL(3,3)=TL(1,2)
      GOTO 37
   34 DO 36 I=2,4
      DO 35 J=1,3
   35 TL(J,I)=0.5
   36 TL(I-1,I)=0.0
   37 CONTINUE
      DMAX=3.0
      DMIN = 2.4
      DO 110 I=1,3
      COSA(I) = COS(CELL(I+3) * D2R)
  110 SINA(I) = SQRT(1.0- COSA(I)**2 )
      COSGS=(COSA(1)*COSA(2)-COSA(3)) / SINA(1)/SINA(2)
      Q(1,1)=CELL(1)*SINA(2)*SQRT(1.0-COSGS**2)
      Q(1,2)=0.0
      Q(1,3)=0.0
      Q(2,1)=-CELL(1)*COSGS*SINA(2)
      Q(2,2)=CELL(2)*SINA(1)
      Q(2,3)=0.0
      Q(3,1)=CELL(1)*COSA(2)
      Q(3,2)=CELL(2)*COSA(1)
      Q(3,3)=CELL(3)
      CALL FILINQ (IATOMS, 'ATOMS', 'FORMATTED', 'INPUT', KINQ)
      IF (KINQ.EQ.-1) CALL KERROR (' No ATOMS file found', 0, 'DISCHK')
      CALL ATOMIN (IATOMS, ATXYZ, ATNAME, IZAT, MAXAT, NA, KEYT)
      NATX = NA
      CALL KERNZA(0.0, XGRAV, 3)
      DO 130 N=1,NA
      DO 130 I=1,3
  130 XGRAV(I) = XGRAV(I) + ATXYZ(I,N)
      DO 131 I=1,3
  131 XGRAV(I) = XGRAV(I) / FLOAT(NA)
      DO 140 IS=1,NPK
      DO 140 I=1,3
      IF (XGRAV(I) + SH(I,IS) .GT. 1.0) SH(I,IS) = SH(I,IS) - 1.0
  140 CONTINUE
      WRITE (LIS2, 149)
  149 FORMAT (/' Shifts, with distance check'/ ' Peak height ', 3X,
     *   'Tx',6X,'Ty',6X,'Tz', 8X, 'contacts by symmetry: Angstrom')
      WRITE (LIS2, 150) DMAX
  150 FORMAT (63X, 'max: ', F4.2)
      K=0
      DO 168 ISS=2,NSYMC
      DO 167 IL=1,NLATT
      DO 166 N=1,NA
      K=K+1
      DO 165 I=1,3
      WORK(I,K)=TL(I,IL)+TT(I,ISS)
      DO 164 J=1,3
  164 WORK(I,K)=WORK(I,K)+ATXYZ(J,N)*IRR(I,J,ISS)
  165 WORK(I,K)=AMOD(4.0+WORK(I,K),1.0)
  166 CONTINUE
  167 CONTINUE
  168 CONTINUE
      RHOMIN = 0.5 * RHO(1)
      NPRINT = 0
      ISAVE = 0
      NBADX = 0
      DO 200 IS=1,NPK
      IBOTS(IS)=0
      BOTS(IS)=999.
      ISAME(IS) = 0
      IF (RHO(IS) .LT. RHOMIN) THEN
         IF (ISAVE .EQ. 0) THEN
            ISAVE = IS - 1
            WRITE (LIS2, FMT='(/
     *        '' The following shift(s) are rejected:'' /)')
            ENDIF
         WRITE (LIS2,2) IS,RHO(IS),(SH(J,IS),J=1,3)
         RHO(IS) = - 1.0
         RHOMIN = RHO(1) + 1.
         GOTO 200
         ENDIF
      RATIO(IS) = 1.
      IF (IS .GT. 1) THEN
         IF (RHO(IS) .GT. RHO(IS-1)) RHO(IS) = RHO(IS-1)
         RATIO(IS) = RHO(IS-1) / RHO(IS)
         IF (RATIO(IS) .GT. 1.25) RHOMIN = RHO(IS) * 0.99
         ENDIF
      RATI2 = 0.0
      IF (IS .GT. 2) RATI2 = RHO(IS-2) / RHO(IS)
      IF (RATI2 .GT. 1.30) RHOMIN = RHO(IS) * 0.99
      RATIO(1) = AMAX1 (RATIO(1), RATIO(IS))
      IF (IS .GE. 5 ) THEN
         TQ = RATIO(IS) + RATIO(IS-1) + RATIO(IS-2) - 3.0
         IF (RATIO(1)-1.0 .GT. 4.0*TQ .AND. RHO(IS) .LT. 0.7*RHO(1)
     *      .AND. IS .LE. 10)  RHOMIN = RHO(IS) * 0.99
         ENDIF
      NPRIN = 0
      DMAX=3.0
      DMAX2 = DMAX * DMAX
      NBAD = 0
      WRITE (LIS2,2) IS,RHO(IS),(SH(J,IS),J=1,3)
    2 FORMAT(I4, F7.0, 3F8.4)
      DO 197 N=1,NA
      DO 170 I=1,3
  170 XY(I)=ATXYZ(I,N)+SH(I,IS)
      IFI=0
      DO 188 ISS=2,NSYMC
      DO 175 I=1,3
      SSH(I)=0.0
      DO 172 J=1,3
  172 SSH(I)=SSH(I)+SH(J,IS)*IRR(I,J,ISS)
  175 SSH(I)=AMOD(4.5+SSH(I),1.0)-0.5
      DO 186 IL=1,NLATT
      IST=IFI+N
      IFI=IFI+NA
      DO 184 J=IST,IFI
      D1=WORK(1,J)+SSH(1)-XY(1)
      D1=(AMOD(4.5+D1,1.0)-0.5)
      D4=ABS(D1)*Q(1,1)
      IF (D4.GT.DMAX) GOTO 184
      D2=WORK(2,J)+SSH(2)-XY(2)
      D2=(AMOD(4.5+D2,1.0)-0.5)
      D5=ABS(D2*Q(2,2)+D1*Q(2,1))
      IF (D5.GT.DMAX) GOTO 184
      D3=WORK(3,J)+SSH(3)-XY(3)
      D3=(AMOD(4.5+D3,1.0)-0.5)
      D6=ABS(D3*Q(3,3)+D2*Q(3,2)+D1*Q(3,1))
      IF (D6.GT.DMAX) GOTO 184
      D = D4*D4+D5*D5+D6*D6
      IF (D.GT.DMAX2) GOTO 184
      D = SQRT(D)
      BOTS(IS) = AMIN1(BOTS(IS), D)
      M=MOD(J-1,NA)+1
      IF ((NPRINT .GT. 50 .AND. NPRIN .EQ. 8) .OR.
     *    (NPRINT .GT. 70 .AND. NPRIN .EQ. 12)) THEN
         DMAX = AMAX1 (BOTS(IS), DMIN)
         DMAX2 = DMAX * DMAX
         WRITE (LIS2, 150) DMAX
         IF (D .GT. DMAX) GOTO 184
         ENDIF
      IF (NPRINT .GT. 99 .AND. NPRIN .EQ. 3) THEN
         DMAX = DMIN
         DMAX2 = DMAX * DMAX
         WRITE (LIS2, 150) DMAX
         IF (D .GT. DMAX) GOTO 184
         ENDIF
      NPRINT = NPRINT + 1
      NPRIN = NPRIN + 1
      IF (NLATT.GT.1) GOTO 180
      WRITE (LIS2,3) ATNAME(N),ATNAME(M),ISS,D
    3 FORMAT(35X, A6, ' ---- ', A6, '  (Sym.', I2, ')  =', F6.2)
      GOTO 182
  180 WRITE (LIS2,4) ATNAME(N),ATNAME(M),ISS,IL,D
    4 FORMAT(35X, A6, ' --- ', A6, ' (Sym.', 2I2, ')  =', F6.2)
  182 CONTINUE
      IF (D .LT. 0.63) ISAME(IS) = ISAME(IS) + 1
      IF (D .GE. 0.43 .AND. D .LE. 0.83) GOTO 193
      IF (N .NE. M) THEN
         IF (MIN0 (IZAT(N), IZAT(M) ) .GE. 7 .AND. D .LE. 2.2) GOTO 193
         IF (D .LE. DMIN) GOTO 193
      ELSE
         WRITE (LIS2, FMT='(63X, ''symmetry?'')')
         ENDIF
  184 CONTINUE
  186 CONTINUE
  188 CONTINUE
      GOTO 197
  193 IBOTS(IS) = IBOTS(IS) + 1
  197 CONTINUE
      IF (IBOTS(IS) .EQ. 0) GOTO 200
      WRITE (LIS2, FMT='(66X, ''? skip'')')
      NBAD = NBAD + 1
      IF (IS .LE. 5) THEN
         NBADX = NBADX + 1
         WRITE (LIS1, 2252) IS, RHO(IS)
 2252    FORMAT (' Warning: shift nr.',I3, ' (pk =', F5.0,
     *   ')  leads to bad contacts')
         IF (ISAME(IS) .GT. 2) WRITE (LIS1, 2253) IS, ISAME(IS)
 2253    FORMAT (' Note for shift nr.',I3, ':', I4,
     *   ' atom pairs coincide by symmetry !')
      ELSE
         IF (NBADX.GT.0) WRITE (LIS1,FMT='(9X,''Also for others ..'')')
         NBADX = -99
         ENDIF
  200 CONTINUE
      IF (ISAVE .GT. 0) NPK = ISAVE
      RETURN
      END
      SUBROUTINE PKSRT (X, MAXAT, IHALF, LIN)
      DIMENSION  X(4, MAXAT)
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     *               WAVE,     CELALL(10),  AMOLW,      ZET,
     *               NELEC,    F000,        ABSMU,      ICENT,
     *               ILATT,    ISYST,       ILAUE,      IMULT,
     *               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     *         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     *         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      COMMON /SEARTA/ D2R, DMPIC, DMAXB, DMOUT, DMINB, ANGM(2), MCON,
     *        SEARDX, NPIC, NATIN, NAT, NATS, NRECY, PSQ, NATREC, KPROG
      PARAMETER (KUSER2=30000)
      COMMON / / NR3D
      INTEGER*2 NR3D(KUSER2)
      INTEGER*2 ILACOM(84000)
      EQUIVALENCE (ILACOM(1), NR3D(1))
      DIMENSION   XS(3), X1(3), IDIFF(19), B(19)
      DIMENSION DXYZM(3)
      DATA MAX, E, D / 0, 0.0, 0.0 /
      DO 101 I = 1, 3
  101 DXYZM(I) = DMPIC * RCELL(I)
      READ (LIN) NNX, NNZ, NNY, NNYT
      NNYOLD=NNY
      IF (IHALF.NE.0) NNY=NNY-3
      NNXP2 = NNX + 2
      NXZ = NNXP2 * (NNZ + 2)
      NXZ3 = 3 * NXZ
      IF (NXZ3 .GT. KUSER2) CALL KERNER (940, 'PKSRT')
      DX = 1.0 / FLOAT(NNX)
      DY = 1.0 / FLOAT(NNYT)
      DZ = 1.0 / FLOAT(NNZ)
      LEVEL = 0
      LIMIT = MIN0(MAXAT, 2*NAT)
 1100 IDIFF(1) = -NXZ - 1
      IDIFF(2) = -NXZ - NNXP2
      IDIFF(3) = -NXZ
      IDIFF(4) = -NXZ + NNXP2
      IDIFF(5) = -NXZ + 1
      IDIFF(6) = -NNXP2 - 1
      IDIFF(7) = -1
      IDIFF(8) = NNXP2 - 1
      IDIFF(9) = -NNXP2
      IDIFF(10) = 0
      DO 1120 I=1,9
      J=20-I
      IDIFF(J) = -IDIFF(I)
 1120 CONTINUE
      NO = 0
      IY = -1
      NY = 0
 1200 REWIND LIN
      READ(LIN)
      READ(LIN)
      IF (IY+2.EQ.NNYOLD) GOTO 1400
      MAX=NXZ
      ISKIP = (NNYOLD-1) * NNZ
      DO 1305 I=1,ISKIP
 1305 READ (LIN)
      CALL RDSECT (MAX, NNXP2, NNZ, NXZ3, LIN)
      REWIND LIN
      READ(LIN)
      READ(LIN)
      CALL RDSECT (MAX, NNXP2, NNZ, NXZ3, LIN)
 1400 MX = MAX - NXZ + NNX + 1
      CALL RDSECT (MAX, NNXP2, NNZ, NXZ3, LIN)
      IY = IY + 1
      NY = MOD(NY+2, 3) - 1
      KK = NXZ3
      IF (NY) 1440, 1460, 1500
 1440 KK = -NXZ3
 1460 DO 1480 I=1,5
      IDIFF(I) = IDIFF(I) - KK
 1480 CONTINUE
      IF (NY .EQ. 0) GO TO 1540
 1500 DO 1520 I=15,19
      IDIFF(I) = IDIFF(I) - KK
 1520 CONTINUE
 1540 DO 2000 IZ=1,NNZ
      MN = MX + 3
      MX = MX + NNXP2
      DO 1980 IX=MN,MX
      IF (NR3D(IX) .LT. LEVEL) GO TO 1980
      DO 1560 I=1,9
      J = IDIFF(I) + IX
      IF (NR3D(IX) .LE. NR3D(J)) GO TO 1980
 1560 CONTINUE
      DO 1580 I=11,19
      J = IDIFF(I) + IX
      IF (NR3D(IX) .LT. NR3D(J)) GO TO 1980
 1580 CONTINUE
      DO 1600 I=1,19
      J = IDIFF(I) + IX
      B(I) = NR3D(J)
 1600 CONTINUE
      B1 = B(3) + B(7) + B(9) + B(11) + B(13) + B(17)
      B2 = B(1) + B(2) + B(4) + B(5) + B(6) + B(8) + B(12) + B(14) +
     +  B(15) + B(16) + B(18) + B(19)
      F = (30.0 * B(10) + 11.0 * B1 - 8.0 * B2) / 63.0
      C = (B(5)+B(12)+B(13)+B(14)+B(19)-B(1)-B(6)-B(7)-B(8)-B(15))/10.0
      DELTAX = C / F
      IF (ABS(DELTAX) .GT. 1.0) GO TO 1620
      D = (B(15)+B(16)+B(17)+B(18)+B(19)-B(1)-B(2)-B(3)-B(4)-B(5))/10.0
      DELTAY = D / F
      IF (ABS(DELTAY) .GT. 1.0) GO TO 1620
      E = (B(4)+B(8)+B(11)+B(14)+B(18)-B(2)-B(6)-B(9)-B(12)-B(16))/10.0
      DELTAZ = E / F
      IF (ABS(DELTAZ) .LE. 1.0) GO TO 1640
 1620 DELTAX = 0.0
      DELTAY = 0.0
      DELTAZ = 0.0
 1640 XX = (FLOAT(IX-MN+1) + DELTAX) * DX
      YY = (FLOAT(IY) + DELTAY) * DY
      ZZ = (FLOAT(IZ) + DELTAZ) * DZ
      A = (9.0 * B(10) + 4.0 * B1 - B2) / 21.0
      BINT = A + 0.5 * (C * DELTAX + D * DELTAY + E * DELTAZ)
      IF (BINT .GT. 1.05 * B(10)) BINT = 1.05 * B(10)
      B(10) =  AMAX1(B(10), BINT)
      NOP1 = NO + 1
      IF(NOP1.GT.MAXAT) GOTO 1821
      X(1,NOP1) = XX
      X(2,NOP1) = YY
      X(3,NOP1) = ZZ
      X(4,NOP1) = B(10)
      IF (NO .EQ. 0) GO TO 1820
      IR=0
      DO 1800 K=1, IMULT
      CALL OPER1 (K, XS, X(1,NOP1))
      DO 1780 I=1,NO
      DO 1720 L=1,3
      X1(L) = X(L,I) - XS(L)
 1680 IF (ABS(X1(L)) .LE. 0.5) GO TO 1700
      X1(L) = X1(L) - SIGN(1.0, X1(L))
      GO TO 1680
 1700 IF (ABS(X1(L)) .GT. DXYZM(L)) GO TO 1780
 1720 CONTINUE
      IF (QUAD2 (X1, X1) .GT. DMPIC) GOTO 1780
      IF (IR.GT.0) X(4,IR)=0.0
      IR=0
      IF (B(10) .LE. X(4,I)) GOTO 1980
      X(1,I) = XX
      X(2,I) = YY
      X(3,I) = ZZ
      X(4,I) = B(10)
      IR=I
 1780 CONTINUE
 1800 CONTINUE
      IF(IR.GT.0) GO TO 1980
 1820 NO = NOP1
 1821 IF (NO .LT. LIMIT) GO TO 1980
      CALL SORT (X, MAXAT, NO, 4)
      NO = NPIC
      LEVEL = X(4,NPIC) + 0.5
 1980 CONTINUE
 2000 CONTINUE
      IF (IY .GE. NNY) GO TO 2100
      IF (IY - NNYOLD + 2) 1400, 1200, 1400
 2100 CALL SORT (X, MAXAT, NO, 4)
      NNN = MIN0 (NO, NPIC)
      IF (NNN .EQ. NPIC) RETURN
      LEVEL = LEVEL - 100
      IF(LEVEL.GE.(-200)) GO TO 1100
      NPIC=NO
      RETURN
      END
      SUBROUTINE PP1
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      EQUIVALENCE (IPR1, IFILE(6)), (LIS2, IFILE(8))
      EQUIVALENCE (IFMAP, IFILE(17)), (ISCRA, IFILE(18))
      EQUIVALENCE (KEYS(27), IMAP),   (KEYS(28), IHALF)
      COMMON /FFTDA/ SCALE, MH(3), NPP(3), XLMIN(3), XLMAX(3)
      DIMENSION ITLE(20)
      EQUIVALENCE (SCALOR,ITLE(18))
      PARAMETER (KUSER2=30000, LUSER2=KUSER2/2)
      COMMON / / X(KUSER2)
      COMPLEX YCOM(LUSER2)
      EQUIVALENCE (X(1), YCOM(1))
      DIMENSION BLACOM(42000)
      EQUIVALENCE (BLACOM(1), X(1))
      INTEGER  P1, P2, R, SKIP, RECS,  D(5)
      NSIZE = KUSER2 - 1
      NX = NPP(1)
      NY = NPP(2)
      NZ = NPP(3)
      RECMAX = 32000.
      MAXSIZ = IFIX ( SQRT (RECMAX * 0.25 * FLOAT (NX*NY*NZ)) )
      NSIZE  = MIN0 (NSIZE, MAXSIZ)
      P1 = NSIZE / (2*NY*(MH(3) + 1))
      P2 = NSIZE / (NX*NZ)
      MPASS = MAX0 (1, NX / MAX0 (1, P1) )
      NBYTES = P1 * P2 * (MH(3)+1) * 8 + 8
      WRITE (LIS2, 101) NBYTES
  101 FORMAT (' Scratch file has max.', I6, ' bytes per record')
      IF (KEYS(26) .EQ. 0 .AND. MPASS .GT. 1) WRITE (IPR1, 102) MPASS
  102 FORMAT (' Appr. number of Fourier-transform passes needed:', I3)
      MPASS = MPASS / MAX0 (1, MPASS/4)
      NPASS = 0
      REWIND ISCRA
      R = -MH(1)
  111 IF (R+P1 .GT. MH(1)) P1 = MH(1) + 1 - R
      CALL RDHKL (YCOM, NY, MH(3)+1, P1, R)
      D(1) = 2*(NY*P1*(MH(3)+1))
      D(2) = 2
      D(3) = D(1)
      D(4) = D(1)
      D(5) = 2*NY
      CALL CMPLFT (X(1), X(2), NSIZE, NY, D)
      CALL WRITEY (YCOM, NY, MH(3)+1, P1, P2, ISCRA)
      NPASS = NPASS + 1
      IF (MPASS .GT. 1 .AND. KEYS(26) .LE. 11 .AND.
     *   MOD (NPASS, MPASS) .EQ. 0) WRITE (IPR1, FMT=
     *   '('' Intermediate transform,  pass'', I2)') NPASS
      R = R + P1
      IF (R .LE. MH(1)) GOTO 111
      IF (IMAP.EQ.5) THEN
         SCALE = 500.0 / SCALE
      ELSE
         SCALE = 3000.0 / SCALE
         ENDIF
      SCALOR = SCALE
      REWIND IFMAP
      NYNEW = NY
      WRITE (IFMAP) ITLE, IMAP, IHALF
      IF (IHALF.NE.0) NYNEW = MIN0 (NY-NY/2+3, NY)
      WRITE (IFMAP) NX,NZ,NYNEW,NY
      REWIND ISCRA
      SKIP = 0
      R = 0
      P1 = NSIZE / (2*NY*(MH(3) + 1))
      RECS = (NY - 1)/P2
      NPASS = 0
  200 IF (R+P2 .GT. NY) P2 = NY - R
      CALL READHL (YCOM, NX,NZ/2, P2,MH(1),MH(3),P1,SKIP,RECS,ISCRA)
      IF (R + P2 .LT. NY) REWIND ISCRA
      SKIP = SKIP + 1
      D(1) = NX * NZ * P2
      D(2) = NZ
      D(3) = NZ * NX
      D(4) = 2 * (MH(3) + 1)
      D(5) = 2
      CALL CMPLFT (X(1), X(2), NSIZE, NX, D)
      D(2) = 2
      D(3) = D(1)
      D(4) = D(1)
      D(5) = NZ
      CALL HERMFT (X(1), X(2), NSIZE, NZ/2, D)
      CALL OUTPUT (X, NZ, NX, P2, R, NY)
      NPASS = NPASS + 1
      IF (MPASS .GT. 1 .AND. KEYS(26) .LE. 11 .AND.
     *   MOD (NPASS, MPASS) .EQ. 0) WRITE (IPR1, FMT=
     *   '('' Final transform, pass'', I2)') NPASS
      R = R + P2
      IF (R .LT. NY) GOTO 200
      CALL FILCLO (ISCRA, 'DELETE')
      RETURN
      END
      SUBROUTINE RDHKL (X, NY, NZ, NX, HS)
      INTEGER HS
      COMPLEX X(NY,NZ,NX)
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      EQUIVALENCE (LIS2, IFILE(8)), (IBINFF, IFILE(16))
      EQUIVALENCE (KEYS(27), IMAP)
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     *               WAVE,     CELALL(10),  AMOLW,      ZET,
     *               NELEC,    F000,        ABSMU,      ICENT,
     *               ILATT,    ISYST,       ILAUE,      IMULT,
     *               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     *         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     *         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      COMMON /FFTDA/ SCALE, MH(3), NPP(3), XLMIN(3), XLMAX(3)
      PARAMETER (MAXBUF = 198)
      DIMENSION FITFFT(5), BUFFFT(MAXBUF)
      EQUIVALENCE (FITFFT(4),EI), (FITFFT(5),PHI)
      INTEGER H, HM, HL
      DIMENSION IHKLX(3,24), PHX(24), NSYMEX(24), GRI(2,24), TAB(15)
      DIMENSION IHKLI(3), ITSYMM(3,24), NRCO(5)
      DATA NEX / -1 /
      DATA NPASS / 0 /
      D2R = ATAN(1.0) / 45.0
      IF (NEX .GE. 0) GOTO 210
      NEX = 0
      SCALE = 0.0
      CALL KERNZI (0, NRCO, 5)
      DO 110 I=1,15
  110 TAB(I) = SIN (FLOAT(30*I) * D2R)
      DO 113 J=1,NSYMM
      DO 113 I=1,3
  113 ITSYMM(I,J) = NINT (TSYMM(I,J) * 12.0)
  210 CALL BINIFF (1, IBINFF, 'BINFFT', FITFFT, NITFFT, BUFFFT, NEND)
      NPASS = NPASS + 1
      HM = HS + NX - 1
      HL = HS
      IDISC = HL * HM
      IHMAX = MAX0 (IABS(HL), IABS(HM))
      IHMIN = MIN0 (IABS(HL), IABS(HM))
      DO 220 H=1,NX
      DO 220 L=1,NZ
      DO 220 K=1,NY
  220 X(K,L,H) = CMPLX(0.0,0.0)
  320 CALL BINIFF (0, IBINFF, 'BINFFT', FITFFT, NITFFT, BUFFFT, NEND)
      IF (NEND.LT.0) GOTO 500
      IF (NPASS .GT. 1) GOTO 330
      NRCO(1) = NRCO(1) + 1
      IF (EI .LT. 0.0) THEN
         NRCO(2) = NRCO(2) + 1
         GOTO 320
         ENDIF
  330 CALL KERF2I(FITFFT, IHKLI, 3)
      IF (ISYST.GT.4 .AND. ISYST.LT.8) GOTO 360
      INMAX = MAX0 (IABS(IHKLI(1)), IABS(IHKLI(2)), IABS(IHKLI(3)))
      INMIN = MIN0 (IABS(IHKLI(1)), IABS(IHKLI(2)), IABS(IHKLI(3)))
      IF (IDISC.LT.0) GOTO 360
      IF (INMIN.GT.IHMAX .OR. INMAX.LT.IHMIN) GOTO 320
  360 CALL FEXPAN (IHKLI, IHKLX, PHX, NSYMEX, NEXP)
      PHI = PHI * D2R
      EC  = EI * COS(PHI)
      ES  = EI * SIN(PHI)
      DO 380 J=1,NEXP
      H = IHKLX(1,J)
      IF (H.GT.HM .OR. H.LT.HL) GOTO 380
      IF (IABS(H) .GE. NPP(1)/2 .OR. IABS(H) .GT. MH(1)) THEN
         NRCO(3) = NRCO(3) + 1
         GOTO 380
         ENDIF
      K = IHKLX(2,J)
      IF (IABS(K) .GE. NPP(2)/2 .OR. IABS(K) .GT. MH(2)) THEN
         NRCO(4) = NRCO(4) + 1
         GOTO 380
         ENDIF
      L = IHKLX(3,J)
      IF (IABS(L) .GE. NPP(3)/2 .OR. IABS(L) .GT. MH(3)) THEN
         NRCO(5) = NRCO(5) + 1
         GOTO 380
         ENDIF
      NU = 0
      IF (IMAP.LT.3.OR.IMAP.GT.4) GOTO 370
      NSYMX = NSYMEX(J)
      NU = - IHKLI(1)*ITSYMM(1,NSYMX) - IHKLI(2)*ITSYMM(2,NSYMX)
     *     - IHKLI(3)*ITSYMM(3,NSYMX)
      NU = MOD(NU,12)
  370 IF (NU.LE.0) NU = NU + 12
      XS = TAB(NU)
      XC = TAB(NU+3)
      GRI(1,J) =  XC*EC - XS*ES
      GRI(2,J) = (XS*EC + XC*ES) * PHX(J)
      NEX = NEX + 1
      SCALE = SCALE + SQRT(GRI(1,J)*GRI(1,J)+GRI(2,J)*GRI(2,J))
      NOKO = 0
      IF (H.EQ.0 .AND. L.EQ.0 .AND. K.NE.0) NOKO = NY - K + 1
      H = H - HS + 1
      IF (K.LT.0) K = NY + K
      K = K + 1
      L = L + 1
      X(K,L,H) = CMPLX(GRI(1,J),GRI(2,J))
      IF (NOKO.NE.0) X(NOKO,L,H) = CONJG(X(K,L,H))
  380 CONTINUE
      GOTO 320
  500 CONTINUE
      IF (NX + HS .LE. MH(1)) RETURN
      CALL FILCLO (IBINFF, 'DELETE')
      WRITE (LIS2, 631) NPASS
  631 FORMAT (' Intermediate transforms required ', I3, ' passes')
      WRITE (LIS2, 690) NRCO(1)
  690 FORMAT (' Number of reflections from input file   =',I7)
      IF (NRCO(2) .GT. 0) WRITE (LIS2, 691) NRCO(2)
  691 FORMAT (' of which',I7,' were rejected'/)
      WRITE (LIS2, 692) NEX
  692 FORMAT (' Number of reflections in one hemisphere =',I7)
      IF (NEX .EQ. 0) CALL KERROR ('No reflections found', 0,'RDHKL')
      IF (NRCO(3).GT.0 .OR. NRCO(4).GT.0 .OR. NRCO(5).GT.0)
     *    WRITE (LIS2, 693)
  693 FORMAT (' not included in calculations, because: '/)
      IF (NRCO(3).GT.0) WRITE (LIS2, 694) MH(1), NRCO(3)
  694 FORMAT (8X,' having H greater than ',I3,'  were ',I6/)
      IF (NRCO(4).GT.0) WRITE (LIS2, 695) MH(2), NRCO(4)
  695 FORMAT (8X,' having K greater than ',I3,'  were ',I6/)
      IF (NRCO(5).GT.0) WRITE (LIS2, 696) MH(3), NRCO(5)
  696 FORMAT (8X,' having L greater than ',I3,'  were ',I6/)
      RETURN
      END
      SUBROUTINE FEXPAN (IHKLI, IHKLX, PHX, NSYMEX, NEXP)
      DIMENSION IHKLI(3), IHKLX(3,24), PHX(24), NSYMEX(24)
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     *               WAVE,     CELALL(10),  AMOLW,      ZET,
     *               NELEC,    F000,        ABSMU,      ICENT,
     *               ILATT,    ISYST,       ILAUE,      IMULT,
     *               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     *         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     *         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      DIMENSION IHKLS(3)
      NEXP = 0
      DO 3900 J=1,NSYMM
      CALL VXMATI (IHKLI, IRSYMM(1,1,J), IHKLS)
      B1 = 1.0
      IF (IHKLS(3)) 1590, 1570, 1600
 1570 IF (IHKLS(1)) 1590, 1580, 1600
 1580 IF (IHKLS(2)) 1590, 1600, 1600
 1590 B1 = -1.0
      DO 3610 I=1,3
 3610 IHKLS(I) = -IHKLS(I)
 1600 CONTINUE
      IF (J.EQ.1) GOTO 1630
      DO 3620 JJ=1, NEXP
      IF (IHKLX(1,JJ) .EQ. IHKLS(1) .AND. IHKLX(2,JJ) .EQ. IHKLS(2)
     *   .AND. IHKLX(3,JJ) .EQ. IHKLS(3)) GOTO 3900
 3620 CONTINUE
 1630 NEXP = NEXP + 1
      CALL KERNAI (IHKLS, IHKLX(1, NEXP), 3)
      PHX(NEXP) = B1
      NSYMEX(NEXP) = J
 3900 CONTINUE
      RETURN
      END
      SUBROUTINE WRITEY (X, NY, NZ, NX, SIZE, ISCRA)
      COMPLEX X(NY,NZ,NX)
      INTEGER SIZE, H, P, Q, R
      P = SIZE
      Q = 0
  100 R = Q + 1
      IF (Q+P .GT. NY) P = NY - Q
      Q = Q + P
      WRITE (ISCRA) (((X(K,L,H), H=1,NX), K=R,Q), L=1,NZ)
      IF (Q .LT. NY) GOTO 100
      RETURN
      END
      SUBROUTINE READHL (X,NX,NZ,NY,HMAX,LMAX,SIZE,SKIP,RECS,ISCRA)
      INTEGER HMAX, SIZE, SKIP ,RECS, H, HL, HU, P, Q
      COMPLEX X(NZ,NX,NY)
      LM = LMAX + 1
      P = SIZE
      HU = NX - HMAX
      IF (SKIP .LE. 0) GOTO 200
      DO 100 Q=1,SKIP
      READ (ISCRA)
  100 CONTINUE
  200 IF (HU + P .GT. NX) GOTO 400
      HL = HU + 1
      HU = HU + P
      READ (ISCRA) (((X(L,H,K), H=HL,HU), K=1,NY), L=1,LM)
      IF (RECS .LE. 0) GOTO 200
      DO 300 Q=1,RECS
      READ (ISCRA)
  300 CONTINUE
      GOTO 200
  400 IF (HU .NE. NX) GOTO 700
      HU = 0
  500 IF (HU+P .GT. HMAX+1) P = HMAX + 1 - HU
      HL = HU + 1
      HU = HU + P
      READ (ISCRA) (((X(L,H,K), H=HL,HU), K=1,NY), L=1,LM)
  550 IF (HU .EQ. HMAX + 1) GOTO 800
      IF (RECS .LE. 0) GOTO 500
      DO 600 Q=1,RECS
      READ (ISCRA)
  600 CONTINUE
      GOTO 500
  700 HL = HU + 1
      HU = HU + P - NX
      IF (HU .GT. HMAX+1) HU = HMAX + 1
      READ (ISCRA) (((X(L,H,K), H=HL,NX), (X(L,H,K), H=1,HU),
     *   K=1,NY), L=1,LM)
      GOTO 550
  800 DO 900 H=2,HU
      HL = NX + 2 - H
      DO 900 K=1,NY
      X(1,HL,K) = CONJG(X(1,H,K))
  900 CONTINUE
      HL = HMAX + 2
      HU = NX - HMAX
      IF (HU .LT. HL) GOTO 920
      DO 910 L=1,LM
      DO 910 K=1,NY
      DO 910 H=HL,HU
      X(L,H,K) = CMPLX(0.0,0.0)
  910 CONTINUE
  920 IF (LM .GE. NZ) GOTO 940
      P = LM + 1
      DO 930 K=1,NY
      DO 930 L=P,NZ
      DO 930 H=1,NX
      X(L,H,K) = CMPLX(0.0,0.0)
  930 CONTINUE
  940 RETURN
      END
      SUBROUTINE OUTPUT (X, NZ, NX, NY, Y, NYT)
      REAL X(NZ,NX,NY)
      INTEGER Y
      COMMON /SYSTA/ IFILE(20), IFSTAT(20), ITIME(4),   KEYS(28),
     *               NFNUM,     NLIT,       NCOLN(32),  NCOLL(32),
     *               NFDOT(32), NFDOL(32),  NLUSER(32), FNUM(32),
     *               SWITCH(28)
      LOGICAL        SWITCH
      COMMON /SYSTB/ PROGNM,    PROSNM,     CCODE,      TITLE,
     *               CHIN,      LIT(32),    CHOUT
      CHARACTER      PROGNM *8, PROSNM *6,  CCODE *6,   TITLE *64,
     *               CHIN *80,  LIT *6,     CHOUT *72
      EQUIVALENCE (LIS2, IFILE(8)), (IFMAP, IFILE(17))
      EQUIVALENCE (KEYS(28), IHALF)
      LOGICAL PRIMAP
      EQUIVALENCE (PRIMAP, SWITCH(11))
      COMMON /FFTDA/ SCALE, MH(3), NPP(3), XLMIN(3), XLMAX(3)
      INTEGER SEC, XL, XU
      INTEGER*2 INUT(250)
      DIMENSION LINE(50)
      DATA NCOL /36/
      IF (IHALF.NE.0)THEN
         NSECO = MIN0 (NYT-NYT/2+2, NYT)
         XLMAX(2) = AMIN1 (XLMAX(2), 0.5)
      ELSE
         NSECO = NYT
         ENDIF
      NXM = MIN0 (NX,  INT(XLMAX(1)*FLOAT(NX)) +1)
      NYM = MIN0 (NYT, INT(XLMAX(2)*FLOAT(NYT))+1)
      NZM = MIN0 (NZ,  INT(XLMAX(3)*FLOAT(NZ)) +1)
      ASC = 1000.0 / FLOAT(NX)
      DO 900 K=1,NY
      SEC = Y + K - 1
      NRSEC = SEC + 1
      DO 350 J=1,NZ
      DO 330 I=1,NX
      X(J,I,K) = X(J,I,K) * SCALE
      INUT4 = MIN0 (NINT (X(J,I,K)), 32767)
  330 INUT(I) = INUT4
      IF (NRSEC.LE.NSECO) WRITE (IFMAP) SEC, J, NX, (INUT(I),I=1,NX)
      IF (NRSEC.EQ.NYT .AND. IHALF.NE.0)
     +WRITE (IFMAP) SEC, J, NX, (INUT(I),I=1,NX)
  350 CONTINUE
      IF (.NOT. PRIMAP) GOTO 900
      XU = 0
  450 XL = XU + 1
      XU = XU + NCOL
      IF (XU.GT.NXM) XU = NXM
      DO 460 L=XL,XU
      INUT(L) = FLOAT((L-1))*ASC + 0.5
  460 CONTINUE
      NEC = 1000.*(FLOAT(SEC))/(FLOAT(NYT)) + 0.5
      WRITE (LIS2, 480) NEC, TITLE
  480 FORMAT ('1SECTION  Y = ', I4, 4X, A64 /)
      WRITE (LIS2, 500) (INUT(NN), NN=XL,XU,2)
  500 FORMAT (8X,'X =',I4,17I6)
      ILIM = XL + 1
      WRITE (LIS2, 520) (INUT(NN), NN=ILIM,XU,2)
  520 FORMAT (12X,18I6)
      WRITE (LIS2, 530)
  530 FORMAT ('0')
      DO 700 I=1,NZM
      LINE(1) = FLOAT((I-1))*1000./FLOAT(NZ) + 0.5
      L = 1
      DO 600 J=XL,XU
      L = L + 1
      LINE(L) = NINT (0.1 * X(I,J,K))
      IF (LINE(L) .LT. -99) LINE(L) = -99
  600 CONTINUE
      WRITE (LIS2, 630) (LINE(J), J=1,L)
  630 FORMAT ('0Z =',I4,' *  ',36I3)
  700 CONTINUE
      IF (XU.LT.NXM)GOTO 450
      IF (SEC.GE.NYM) PRIMAP = .FALSE.
  900 CONTINUE
      RETURN
      END
      SUBROUTINE CMPLFT (X, Y, NSIZE, N, D)
      REAL X(NSIZE), Y(NSIZE)
      INTEGER D(5),PMAX,PSYM,TWOGRP,FACTOR(15),SYM(15),UNSYM(15)
      PMAX   = 5
      TWOGRP = 4
      CALL SRFP (N, PMAX, TWOGRP, FACTOR, SYM, PSYM, UNSYM)
      CALL MDFTKD (N, FACTOR, D, X, Y, NSIZE)
      CALL DIPRP (N, SYM, PSYM, UNSYM, D, X, Y, NSIZE)
      RETURN
      END
      SUBROUTINE SRFP (PTS,PMAX,TWOGRP,FACTOR,SYM,PSYM,UNSYM)
      INTEGER PTS,PMAX,TWOGRP,PSYM, FACTOR(15), SYM(15), UNSYM(15)
      INTEGER PP(14), QQ(7), F, P, PTWO, Q, R
      N = PTS
      PSYM = 1
      F = 2
      P = 0
      Q = 0
  100 IF (N.LE.1) GOTO 500
      DO 200 J=F,PMAX
      IF (N.EQ.(N/J)*J) GOTO 300
  200 CONTINUE
      CALL KERNER (200, 'SRFP')
  300 F = J
      N = N / F
      IF (N.EQ.(N/F)*F) GOTO 400
      Q = Q + 1
      QQ(Q) = F
      GOTO 100
  400 N = N / F
      P = P + 1
      PP(P) = F
      PSYM = PSYM * F
      GOTO 100
  500 R = 1
      IF (Q.EQ.0) R = 0
      IF (P.LT.1) GOTO 700
      DO 600 J=1,P
      JJ = P + 1 - J
      SYM(J) = PP(JJ)
      FACTOR(J) = PP(JJ)
      JJ = P + Q + J
      FACTOR(JJ) = PP(J)
      JJ = P + R + J
      SYM(JJ) = PP(J)
  600 CONTINUE
  700 IF (Q.LT.1) GOTO 900
      DO 800 J=1,Q
      JJ = P + J
      UNSYM(J) = QQ(J)
      FACTOR(JJ) = QQ(J)
  800 CONTINUE
      SYM(P+1) = PTS / PSYM**2
  900 JJ = 2*P + Q
      FACTOR(JJ+1) = 0
      PTWO = 1
      J = 0
 1000 J = J + 1
      IF (FACTOR(J).EQ.0) GOTO 1200
      IF (FACTOR(J).NE.2) GOTO 1000
      PTWO = PTWO * 2
      FACTOR(J) = 1
      IF (PTWO .GE. TWOGRP) GOTO 1100
      IF (FACTOR(J+1).EQ.2) GOTO 1000
 1100 FACTOR(J) = PTWO
      PTWO = 1
      GOTO 1000
 1200 IF (P.EQ.0) R = 0
      JJ = 2*P + R
      SYM(JJ+1) = 0
      IF (Q.LE.1) Q = 0
      UNSYM(Q+1) = 0
      RETURN
      END
      SUBROUTINE DIPRP (PTS, SYM, PSYM, UNSYM, DIM, X, Y, NSIZE)
      REAL X(NSIZE), Y(NSIZE)
      INTEGER SYM(15), UNSYM(15), DIM(5), PTS, PSYM, DK, PUNSYM, TEST
      LOGICAL ONEMOD
      INTEGER SEP, DELTA, P, P0, P1, P2, P3, P4, P5, SIZE
      INTEGER V(14), MODULO(14), S(14), U(14)
      DATA MODS / 0 /
      NEST = 14
      NT = DIM(1)
      SEP = DIM(2)
      P2 = DIM(3)
      SIZE = DIM(4) - 1
      P4 = DIM(5)
      IF (SYM(1).EQ.0) GOTO 500
      DO 100 J=1,NEST
      U(J) = 1
      S(J) = 1
  100 CONTINUE
      N = PTS
      DO 200 J=1,NEST
      IF (SYM(J).EQ.0) GOTO 300
      JJ = NEST + 1 - J
      U(JJ) = N
      N = N / SYM(J)
      S(JJ) = N
  200 CONTINUE
  300 JJ = 0
      L = 1
      V(1) = 1
  310 L = L + 1
      V(L) = V(L-1)
  320 IF (L.LT.NEST) GOTO 310
      N = V(NEST)
      JJ = JJ + 1
      IF (JJ.GE.N) GOTO 400
      DELTA = (N-JJ) * SEP
      P1 = (JJ-1)*SEP + 1
      DO 350 P0=P1,NT,P2
      P3 = P0 + SIZE
      DO 350 P=P0,P3,P4
      P5 = P + DELTA
      T = X(P)
      X(P) = X(P5)
      X(P5) = T
      T = Y(P)
      Y(P) = Y(P5)
      Y(P5) = T
  350 CONTINUE
  400 V(L) = V(L) + S(L)
      IF (V(L).LE.U(L)) GOTO 320
      L = L - 1
      IF (L.NE.0) GOTO 400
  500 IF (UNSYM(1).EQ.0) GOTO 1900
      PUNSYM = PTS / PSYM**2
      MULT = PUNSYM / UNSYM(1)
      TEST = (UNSYM(1)*UNSYM(2)-1) * MULT * PSYM
      LK = MULT
      DK = MULT
      DO 600 K=2,NEST
      IF (UNSYM(K).EQ.0) GOTO 700
      LK = LK * UNSYM(K-1)
      DK = DK / UNSYM(K)
      U(K) = (LK-DK) * PSYM
      MODS = K
  600 CONTINUE
  700 ONEMOD = MODS.LT.3
      IF (ONEMOD) GOTO 900
      DO 800 J=3,MODS
      JJ = MODS + 3 - J
      MODULO(JJ) = U(J)
  800 CONTINUE
  900 MODULO(2) = U(2)
      JL = (PUNSYM-3) * PSYM
      MS = PUNSYM * PSYM
      DO 1800 J=PSYM,JL,PSYM
      K = J
 1000 K = K * MULT
      IF (ONEMOD) GOTO 1200
      DO 1100 I=3,MODS
      K = K - (K/MODULO(I))*MODULO(I)
 1100 CONTINUE
 1200 IF (K.GE.TEST) GOTO 1300
      K = K - (K/MODULO(2))*MODULO(2)
      GOTO 1400
 1300 K = K - (K/MODULO(2))*MODULO(2)+MODULO(2)
 1400 IF (K.LT.J) GOTO 1000
      IF (K.EQ.J) GOTO 1800
      DELTA = (K-J) * SEP
      DO 1600 L=1,PSYM
      DO 1500 M=L,PTS,MS
      P1 = (M+J-1)*SEP + 1
      DO 1500 P0=P1,NT,P2
      P3 = P0 + SIZE
      DO 1500 JJ=P0,P3,P4
      KK = JJ + DELTA
      T = X(JJ)
      X(JJ) = X(KK)
      X(KK) = T
      T = Y(JJ)
      Y(JJ) = Y(KK)
      Y(KK) = T
 1500 CONTINUE
 1600 CONTINUE
 1800 CONTINUE
 1900 RETURN
      END
      SUBROUTINE MDFTKD (N, FACTOR, DIM, X, Y, NSIZE)
      INTEGER FACTOR(15), DIM(5), F, P, R, S
      REAL X(NSIZE), Y(NSIZE)
      S = DIM(2)
      F = 0
      M = N
  100 F = F + 1
      P = FACTOR(F)
      IF (P.EQ.0) RETURN
      M = M / P
      R = M * S
      NDIR1 = NSIZE - R
      NDIR2 = NDIR1 - R
      NDIR3 = NDIR2 - R
      NDIR4 = NDIR3 - R
      GOTO (100, 200, 300, 400, 500), P
  200 CALL R2CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), DIM, NSIZE, NDIR1)
      GOTO 100
  300 CONTINUE
      CALL R3CFTK(N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1)
     ., DIM, NSIZE, NDIR1, NDIR2)
      GOTO 100
  400 CALL R4CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1)
     ., X(3*R+1), Y(3*R+1), DIM, NSIZE, NDIR1, NDIR2, NDIR3)
      GOTO 100
  500 CALL R5CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1)
     ., X(3*R+1), Y(3*R+1), X(4*R+1), Y(4*R+1), DIM, NSIZE, NDIR1,
     +  NDIR2, NDIR3, NDIR4)
      GOTO 100
      END
      SUBROUTINE R2CFTK (N, M, X0, Y0, X1, Y1, DIM, NDIR0, NDIR1)
      INTEGER DIM(5), SIZE, SEP
      REAL X0(NDIR0), Y0(NDIR0), X1(NDIR1), Y1(NDIR1), IS, IU
      LOGICAL FOLD,ZERO
      DATA TWOPI / 6.2831853 /
      DATA C, S / 0.0, 0.0 /
      NT = DIM(1)
      SEP = DIM(2)
      L1 = DIM(3)
      SIZE = DIM(4) - 1
      K2 = DIM(5)
      NS = N * SEP
      M2 = M * 2
      FM2 = FLOAT(M2)
      MOVER2 = M/2 + 1
      MM2 = SEP * M2
      FJM1 = -1.0
      DO 600 J=1, MOVER2
      FOLD = J.GT.1 .AND. 2*J.LT.M+2
      K0 = (J-1)*SEP + 1
      FJM1 = FJM1 + 1.0
      ANGLE = TWOPI * FJM1 / FM2
      ZERO = ANGLE.EQ.0.0
      IF (ZERO) GOTO 200
      C = COS(ANGLE)
      S = SIN(ANGLE)
      GOTO 200
  100 FOLD = .FALSE.
      K0 = (M+1-J)*SEP + 1
      C = -C
  200 DO 500 KK=K0,NS,MM2
      DO 440 L=KK,NT,L1
      K1 = L + SIZE
      DO 420 K=L,K1,K2
      RS = X0(K) + X1(K)
      IS = Y0(K) + Y1(K)
      RU = X0(K) - X1(K)
      IU = Y0(K) - Y1(K)
      X0(K) = RS
      Y0(K) = IS
      IF (ZERO) GOTO 300
      X1(K) = RU*C + IU*S
      Y1(K) = IU*C - RU*S
      GOTO 420
  300 X1(K) = RU
      Y1(K) = IU
  420 CONTINUE
  440 CONTINUE
  500 CONTINUE
      IF (FOLD) GOTO 100
  600 CONTINUE
      RETURN
      END
      SUBROUTINE R3CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, DIM,
     +     NDIR0, NDIR1, NDIR2)
      REAL X0(NDIR0),Y0(NDIR0),X1(NDIR1),Y1(NDIR1),X2(NDIR2),Y2(NDIR2),
     +      I0,I1,I2,IA,IB,IS
      LOGICAL FOLD,ZERO
      INTEGER DIM(5), SIZE, SEP
      DATA TWOPI / 6.2831853 /
      DATA A/-0.5/, B/0.866/
      DATA C1, S1, C2, S2 / 0.0, 0.0, 0.0, 0.0 /
      NT = DIM(1)
      SEP = DIM(2)
      L1 = DIM(3)
      SIZE = DIM(4) - 1
      K2 = DIM(5)
      NS = N * SEP
      M3 = M * 3
      FM3 = FLOAT(M3)
      MM3 = SEP * M3
      MOVER2 = M/2 + 1
      FJM1 = -1.0
      DO 600 J=1, MOVER2
      FOLD = J.GT.1 .AND. 2*J.LT.M+2
      K0 = (J-1)*SEP + 1
      FJM1 = FJM1 + 1.0
      ANGLE = TWOPI * FJM1 / FM3
      ZERO = ANGLE.EQ.0.0
      IF (ZERO) GOTO 200
      C1 = COS(ANGLE)
      S1 = SIN(ANGLE)
      C2 = C1*C1 - S1*S1
      S2 = S1*C1 + C1*S1
      GOTO 200
  100 FOLD = .FALSE.
      K0 = (M+1-J)*SEP + 1
      T = C1*A + S1*B
      S1 = C1*B - S1*A
      C1 = T
      T = C2*A - S2*B
      S2 = -C2*B - S2*A
      C2 = T
  200 DO 500 KK=K0,NS,MM3
      DO 440 L=KK,NT,L1
      K1 = L + SIZE
      DO 420 K=L,K1,K2
      R0 = X0(K)
      I0 = Y0(K)
      RS = X1(K) + X2(K)
      IS = Y1(K) + Y2(K)
      X0(K) = R0 + RS
      Y0(K) = I0 + IS
      RA = R0 + RS*A
      IA = I0 + IS*A
      RB = (X1(K)-X2(K)) * B
      IB = (Y1(K)-Y2(K)) * B
      IF (ZERO) GOTO 300
      R1 = RA + IB
      I1 = IA - RB
      R2 = RA - IB
      I2 = IA + RB
      X1(K) = R1*C1 + I1*S1
      Y1(K) = I1*C1 - R1*S1
      X2(K) = R2*C2 + I2*S2
      Y2(K) = I2*C2 - R2*S2
      GOTO 420
  300 X1(K) = RA + IB
      Y1(K) = IA - RB
      X2(K) = RA - IB
      Y2(K) = IA + RB
  420 CONTINUE
  440 CONTINUE
  500 CONTINUE
      IF (FOLD) GOTO 100
  600 CONTINUE
      RETURN
      END
      SUBROUTINE R4CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, X3, Y3, DIM,
     +   NDIR0, NDIR1, NDIR2, NDIR3)
      REAL X0(NDIR0),Y0(NDIR0),X1(NDIR1),Y1(NDIR1),
     +     X2(NDIR2),Y2(NDIR2),X3(NDIR3),Y3(NDIR3)
      INTEGER DIM(5), SIZE, SEP
      LOGICAL FOLD,ZERO
      REAL I1,I2,I3,IS0,IS1,IU0,IU1
      DATA TWOPI / 6.2831853 /
      DATA C1, S1, C2, S2, C3, S3 / 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /
      NT = DIM(1)
      SEP = DIM(2)
      L1 = DIM(3)
      SIZE = DIM(4) - 1
      K2 = DIM(5)
      NS = N * SEP
      M4 = M * 4
      FM4 = FLOAT(M4)
      MM4 = SEP * M4
      MOVER2 = M/2 + 1
      FJM1 = -1.0
      DO 600 J=1, MOVER2
      FOLD = J.GT.1 .AND. 2*J.LT.M+2
      K0 = (J-1)*SEP + 1
      FJM1 = FJM1 + 1.0
      ANGLE = TWOPI * FJM1 / FM4
      ZERO = ANGLE.EQ.0.0
      IF (ZERO) GOTO 200
      C1 = COS(ANGLE)
      S1 = SIN(ANGLE)
      C2 = C1*C1 - S1*S1
      S2 = S1*C1 + C1*S1
      C3 = C2*C1 - S2*S1
      S3 = S2*C1 + C2*S1
      GOTO 200
  100 FOLD = .FALSE.
      K0 = (M+1-J)*SEP + 1
      T = C1
      C1 = S1
      S1 = T
      C2 = -C2
      T = C3
      C3 = -S3
      S3 = -T
  200 DO 500 KK=K0,NS,MM4
      DO 440 L=KK,NT,L1
      K1 = L + SIZE
      DO 420 K=L,K1,K2
      RS0 = X0(K) + X2(K)
      IS0 = Y0(K) + Y2(K)
      RU0 = X0(K) - X2(K)
      IU0 = Y0(K) - Y2(K)
      RS1 = X1(K) + X3(K)
      IS1 = Y1(K) + Y3(K)
      RU1 = X1(K) - X3(K)
      IU1 = Y1(K) - Y3(K)
      X0(K) = RS0 + RS1
      Y0(K) = IS0 + IS1
      IF (ZERO) GOTO 300
      R1 = RU0 + IU1
      I1 = IU0 - RU1
      R2 = RS0 - RS1
      I2 = IS0 - IS1
      R3 = RU0 - IU1
      I3 = IU0 + RU1
      X2(K) = R1*C1 + I1*S1
      Y2(K) = I1*C1 - R1*S1
      X1(K) = R2*C2 + I2*S2
      Y1(K) = I2*C2 - R2*S2
      X3(K) = R3*C3 + I3*S3
      Y3(K) = I3*C3 - R3*S3
      GOTO 420
  300 X2(K) = RU0 + IU1
      Y2(K) = IU0 - RU1
      X1(K) = RS0 - RS1
      Y1(K) = IS0 - IS1
      X3(K) = RU0 - IU1
      Y3(K) = IU0 + RU1
  420 CONTINUE
  440 CONTINUE
  500 CONTINUE
      IF (FOLD) GOTO 100
  600 CONTINUE
      RETURN
      END
      SUBROUTINE R5CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, X3, Y3, X4, Y4,
     *                   DIM, NDIR0, NDIR1, NDIR2, NDIR3, NDIR4)
      REAL X0(NDIR0),Y0(NDIR0),X1(NDIR1),Y1(NDIR1),X2(NDIR2),
     +     Y2(NDIR2),X3(NDIR3),Y3(NDIR3),X4(NDIR4),Y4(NDIR4),
     +     I0,I1,I2,I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2
      INTEGER DIM(5), SIZE, SEP
      LOGICAL FOLD,ZERO
      DATA TWOPI / 6.2831853 /
      DATA A1/0.30902/   ,B1/0.95106/   ,A2/-0.80902/   ,B2/0.58779/
      DATA C1, S1, C2, S2, C3, S3, C4, S4 / 0.,0.,0.,0.,0.,0.,0.,0. /
      NT = DIM(1)
      SEP = DIM(2)
      L1 = DIM(3)
      SIZE = DIM(4) - 1
      K2 = DIM(5)
      NS = N * SEP
      M5 = M * 5
      FM5 = FLOAT(M5)
      MM5 = SEP*M5
      MOVER2 = M/2 + 1
      FJM1 = -1.0
      DO 600 J=1, MOVER2
      FOLD = J.GT.1 .AND. 2*J.LT.M+2
      K0 = (J-1)*SEP + 1
      FJM1 = FJM1 + 1.0
      ANGLE = TWOPI * FJM1 / FM5
      ZERO = ANGLE.EQ.0.0
      IF (ZERO) GOTO 200
      C1 = COS(ANGLE)
      S1 = SIN(ANGLE)
      C2 = C1*C1 - S1*S1
      S2 = S1*C1 + C1*S1
      C3 = C2*C1 - S2*S1
      S3 = S2*C1 + C2*S1
      C4 = C2*C2 - S2*S2
      S4 = S2*C2 + C2*S2
      GOTO 200
  100 FOLD = .FALSE.
      K0 = (M+1-J)*SEP + 1
      T = C1*A1 + S1*B1
      S1 = C1*B1 - S1*A1
      C1 = T
      T = C2*A2 + S2*B2
      S2 = C2*B2 - S2*A2
      C2 = T
      T = C3*A2 - S3*B2
      S3 = -C3*B2 - S3*A2
      C3 = T
      T = C4*A1 - S4*B1
      S4 = -C4*B1 - S4*A1
      C4 = T
  200 DO 500 KK=K0,NS,MM5
      DO 440 L=KK,NT,L1
      K1 = L + SIZE
      DO 420 K=L,K1,K2
      R0 = X0(K)
      I0 = Y0(K)
      RS1 = X1(K) + X4(K)
      IS1 = Y1(K) + Y4(K)
      RU1 = X1(K) - X4(K)
      IU1 = Y1(K) - Y4(K)
      RS2 = X2(K) + X3(K)
      IS2 = Y2(K) + Y3(K)
      RU2 = X2(K) - X3(K)
      IU2 = Y2(K) - Y3(K)
      X0(K) = R0 + RS1+RS2
      Y0(K) = I0 + IS1+IS2
      RA1 = R0 + RS1*A1+RS2*A2
      IA1 = I0 + IS1*A1+IS2*A2
      RA2 = R0 + RS1*A2+RS2*A1
      IA2 = I0 + IS1*A2+IS2*A1
      RB1 = RU1*B1 + RU2*B2
      IB1 = IU1*B1 + IU2*B2
      RB2 = RU1*B2 - RU2*B1
      IB2 = IU1*B2 - IU2*B1
      IF (ZERO) GOTO 300
      R1 = RA1 + IB1
      I1 = IA1 - RB1
      R2 = RA2 + IB2
      I2 = IA2 - RB2
      R3 = RA2 - IB2
      I3 = IA2 + RB2
      R4 = RA1 - IB1
      I4 = IA1 + RB1
      X1(K) = R1*C1 + I1*S1
      Y1(K) = I1*C1 - R1*S1
      X2(K) = R2*C2 + I2*S2
      Y2(K) = I2*C2 - R2*S2
      X3(K) = R3*C3 + I3*S3
      Y3(K) = I3*C3 - R3*S3
      X4(K) = R4*C4 + I4*S4
      Y4(K) = I4*C4 - R4*S4
      GOTO 420
  300 X1(K) = RA1 + IB1
      Y1(K) = IA1 - RB1
      X2(K) = RA2 + IB2
      Y2(K) = IA2 - RB2
      X3(K) = RA2 - IB2
      Y3(K) = IA2 + RB2
      X4(K) = RA1 - IB1
      Y4(K) = IA1 + RB1
  420 CONTINUE
  440 CONTINUE
  500 CONTINUE
      IF (FOLD) GOTO 100
  600 CONTINUE
      RETURN
      END
      SUBROUTINE HERMFT (X, Y, NSIZE, N, DIM)
      REAL X(NSIZE), Y(NSIZE)
      INTEGER DIM(5), D2, D3, D4, D5
      DATA TWOPI / 6.2831853 /
      TWON = FLOAT(2*N)
      NT = DIM(1)
      D2 = DIM(2)
      D3 = DIM(3)
      D4 = DIM(4) - 1
      D5 = DIM(5)
      DO 100 I0=1,NT,D3
      I1 = I0 + D4
      DO 100 I=I0,I1,D5
      A = X(I)
      B = Y(I)
      X(I) = A + B
      Y(I) = A - B
  100 CONTINUE
      NOVER2 = N/2 + 1
      IF (NOVER2 .LT. 2) GOTO 500
      DO 400 I0 = 2, NOVER2
      ANGLE = TWOPI * FLOAT(I0-1) / TWON
      CO = COS(ANGLE)
      SI = SIN(ANGLE)
      K = (N + 2 - 2*I0)*D2
      K1 = (I0 - 1)*D2 + 1
      DO 300 I1=K1,NT,D3
      I2 = I1 + D4
      DO 200 I=I1,I2,D5
      J = I + K
      A = X(I) + X(J)
      B = X(I) - X(J)
      C = Y(I) + Y(J)
      D = Y(I) - Y(J)
      E = B*CO + C*SI
      F = B*SI - C*CO
      X(I) = A + F
      X(J) = A - F
      Y(I) = E + D
      Y(J) = E - D
  200 CONTINUE
  300 CONTINUE
  400 CONTINUE
      CALL CMPLFT (X, Y, NSIZE, N, DIM)
  500 RETURN
      END
      SUBROUTINE RDSECT (MAX, NNXP2, NNZ, NXZ3, LIN)
      PARAMETER (KUSER2=30000)
      COMMON / / NR3D
      INTEGER*2 NR3D(KUSER2)
      INTEGER*2 ILACOM(84000)
      EQUIVALENCE (ILACOM(1), NR3D(1))
      IF (MAX .GE. NXZ3) MAX = 0
      MX = MAX
      MAX = MAX - 2
      DO 1320 IZ=1,NNZ
      MIN = MAX + 3
      MAX = MAX + NNXP2
      READ (LIN) IYSEC, IZLIN, IXTOT, (NR3D(IX),IX=MIN,MAX)
      NR3D(MAX+1) = NR3D(MIN)
      NR3D(MAX+2) = NR3D(MIN+1)
 1320 CONTINUE
      MIN = MAX + 3
      MAX = MAX + NNXP2 + NNXP2 + 2
      DO 1340 IX=MIN,MAX
      MX = MX + 1
      NR3D(IX) = NR3D(MX)
 1340 CONTINUE
      RETURN
      END
      SUBROUTINE SORT (X, MAXAT, NAT, N)
      DIMENSION X(4,MAXAT), T(4)
      INT=2
 1000 INT=INT+INT
      IF(INT.LT.NAT)GO TO 1000
      INT = MIN0 (NAT, (3*INT)/4-1)
 1020 INT=INT/2
      IFIN=NAT-INT
      DO 1200 II=1,IFIN
      I=II
      J=I+INT
      IF (X(N,I) .GE. X(N,J)) GOTO 1200
      DO 1060 K=1,4
      T(K) = X(K,J)
 1060 CONTINUE
 1080 DO 1100 K=1,4
      X(K,J) = X(K,I)
 1100 CONTINUE
      J=I
      I=I-INT
      IF (I) 1140, 1140, 1120
 1120 IF (X(N,I) .LT. T(N)) GOTO 1080
 1140 DO 1160 K=1,4
      X(K,J) = T(K)
 1160 CONTINUE
 1200 CONTINUE
      IF(INT.NE.1)GO TO 1020
      RETURN
      END
      FUNCTION QUAD2 (X1, X2)
      DIMENSION X1(3), X2(3)
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     +               WAVE,     CELALL(10),  AMOLW,      ZET,
     +               NELEC,    F000,        ABSMU,      ICENT,
     +               ILATT,    ISYST,       ILAUE,      IMULT,
     +               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     +         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     +         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      CALL VMATV1 (X1, RRMAT, X2, DIST2)
      QUAD2 = DIST2
      RETURN
      END
      SUBROUTINE OPER1 (J, XN, XYZ)
      DIMENSION XN(3), XYZ(3)
      DIMENSION FS(3,3,24)
      COMMON /CRYSA/ CELL(6),  CELLSD(6),   RCELL(6),   VOLUM,
     *               WAVE,     CELALL(10),  AMOLW,      ZET,
     *               NELEC,    F000,        ABSMU,      ICENT,
     *               ILATT,    ISYST,       ILAUE,      IMULT,
     *               IUNIQ,    IPOLA,       NTYPE,      NSYMM,
     *         IRSYMM(3,3,24), TSYMM(3,24), NLATT,      TLATT(3,4),
     *         FRAC2C(3,3),    CART2F(3,3), RRMAT(3,3), SSMAT(3,3)
      LOGICAL CONT
      DATA CONT / .FALSE. /
      IF (CONT) GOTO 111
      CONT = .TRUE.
      CALL KERI2F (IRSYMM, FS, 9*NSYMM)
  111 ISYM = MOD(J,NSYMM)
      IF (ISYM .EQ. 0) ISYM = NSYMM
      IP = (J-1) / NSYMM
      ILAT = MOD(IP,NLATT) + 1
      IF (IP .LT. NLATT) THEN
         DO 120 I = 1,3
  120    XN(I) = XYZ(1) * FS(I,1,ISYM)
     *         + XYZ(2) * FS(I,2,ISYM)
     *         + XYZ(3) * FS(I,3,ISYM) + TSYMM(I,ISYM) + TLATT(I,ILAT)
      ELSE
         DO 130 I = 1,3
  130    XN(I) =-XYZ(1) * FS(I,1,ISYM)
     *         - XYZ(2) * FS(I,2,ISYM)
     *         - XYZ(3) * FS(I,3,ISYM) - TSYMM(I,ISYM) + TLATT(I,ILAT)
         ENDIF
      RETURN
      END
