C
C BK ACHTUNG! REAL*8 und CHARACTER*8 muessen intern 8 Bytes belegen!
C BK          REAL   und CHARACTER*4 muessen intern 4 Bytes belegen!
C
      program dls76
C          DATA SET DHOFCSLS76 AT LEVEL 004 AS OF 12/22/77
C          DATA SET DHOFDLS76  AT LEVEL 001 AS OF 04/06/77
C     PROGRAM DLS76 (INPUT,OUTPUT,PUNCH,TAPE5=INPUT,TAPE6=OUTPUT,TAPE7=P
C    1UNCH,TAPE8,TAPE9,TAPE10)
C
C                   DISTANCE LEAST SQUARES PROGRAM DLS-76
C                   -------------------------------------
C
C                                    BY
C
C                         CH BAERLOCHER AND A HEPP
C                   INSTITUT FUER KRISTALLOGRAPHIE DER ETH
C                        CH-8092 ZURICH, SWITZERLAND
C
C        THIS PROGRAM IS A REVISED AND EXTENDED VERSION OF THE DLS-74
C                                    BY
C                               BRUNO GUIGAS
C                INSTITUT FUER KRISTALLOGRAPHIE DER UNIVERSITAET
C                        D 7500 KARLSRUHE, GERMANY
C
C     IBM VERSION BY S.L. LAWTON AND E.L. WU
C     MOBIL RESEARCH AND DEVELOPMENT CORP., PAULSBORO, N.J. USA
C
C     MODIFIED 1993-04-07 BY BERNT KARASCH FOR TRANSLATION BY F2C.
C     MODIFICATIONS MARKED WITH 'BK'.
C     INSTITUT FUER MINERALOGIE, RUHR-UNIVERSITAET BOCHUM, GERMANY
C     EMAIL: HERMANN.GIES@RUBA.RZ.RUHR-UNI-BOCHUM.DBP.DE.
C
C
C **************************************************************************
C					PORTED TO COMPAQ FORTRAN FOR PC 
C				RJ ANGEL VT CRYSTALLOGRAPHY LAB JAN 2002
C
C  NOTES: THIS VERSION I CALLED DLS76-02, ONLY HAS CHANGES FOR PORTING
C		LATER VERSION(S) WITH LETTERS IN NAME HAVE FURTHER DEVELOPMENTS
C
C	1: INTRODUCED NEW DATE ROUTINE
C	2: NEW SUBROUTINE OPENFILES TO GET NAMES OF FILES AND OPEN THEM
C	3: CHANGED MANY OUTPUT FORMATS (NO CC)
C
C	
C **************************************************************************
      CHARACTER*2 ATID
      CHARACTER*72 ANS,TIT1,TIT2,TIT3,TIT4,P10FN
      INTEGER BT
      DIMENSION IXYZ(3),A1(3)
      REAL*8 ID
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK4/ RN,RSIGMA,RWERT,ICVTST
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200) ,XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK8/ SI(3,150),SIGN(3,3,150),SYMOPS(10,200),B(3,3,200),
     1BI(3,200),IND,MTR(200)
      COMMON /BLOCK9/ DELV(300),IFELD(6),IV(300),IV1(300),KVAR(3,200)
      COMMON S(50000)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
      COMMON /BLOC13/ ISYS,CVGTST,LSYS
      COMMON /BLOC15/ NAPID,IRNG,KRNG
      COMMON /BLOC16/ ITYPE(200),BT(15),IA,NZA,ICON(40,11),MD(40,18),NCA
     1PID,CAPID(15,4),IAT,IAB,IBT(40,4),DOT,WEIGHT(15,3),DELTAX(3),DELTA
     20(3),DA,DB,ANG,TOT(40,4),OTO(40,6),PI180,IMOV(3),MSYMAB(200),DOBIN
     3(400),KTYPE(20),IDK
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
C BK     1IDIMM
      COMMON /BLOC18/ ROV,NC,NA,NS,IVF,IRNG1
C BK
      COMMON /BLOCKBK/ NUMB(10)
      DATA IXYZ/'X','Y','Z'/
C BK B
      DATA IBKSI /'SI'/
      DATA IBKAL /'AL'/
      DATA IBKO  /'O'/
      DATA IBKNA /'NA'/
      DATA IBKP  /'P'/
      DATA NUMB/1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H0/
C BK E
      KRNG=0
      IRNG=0
      IRNG1=0
      NTOUT=12
      KRONEK(1,1)=1
      KRONEK(1,2)=0
      KRONEK(2,1)=0
      KRONEK(2,2)=1
      KRONEK(1,3)=0
      KRONEK(2,3)=0
      KRONEK(3,1)=0
      KRONEK(3,2)=0
      KRONEK(3,3)=1
      KRON(1,1)=1
      KRON(1,2)=-1
      KRON(2,1)=-1
      KRON(2,2)=1
C     CALL TIMER(ITIME1)
C     TIME=ITIME1/38450
      NTIME=0
      IADAPT=0
      IDOB1=0
      PI180=3.1415926535898/180.
C
C  NEW SUBROUTINE CALL TO OPEN FILES AND PRINT HEADER RJA JAN 02
C
	CALL OPENFILES
C
    1 CONTINUE
      NCAPID=0
      RWM=0.5
      ICVTST=0
      NCYC=0
C***  PRINT EXECUTION TIME (CPU).
C       GIVES THE ELAPSED TOTAL TIME.
C     ITIME2=ITIME1
C     CALL TIMER(ITIME1)
C     TIME=ITIME1/38450
C     IF (NTIME.EQ.0) GO TO 111
C     TIME=(ITIME1-ITIME2)/38450
C     WRITE (NTOUT,1110) TIME
C1110 FORMAT (16H0EXECUTION TIME:,3H = ,F15.3,4H SEC)
C*** END OF TIME PRINTING
C
C      WRITE (NTOUT,1008)
      IZEIL=IZEIL+2
C     CALL KOPF1
C 1008 FORMAT (1X,12(6H------))
  111 NTIME=1
   40 CALL DATIN
   41 IF ((IRNG.GT.0).AND.(KRNG.GT.1)) GO TO 12
      CALL DATOUT
      IF (IERR.GT.0) GO TO 1
      GO TO 10
   12 IF (IZEIL.GT.0) CALL KOPF
      CALL PAROUT
   10 NNCY=NCYCLE+1
      NCB=1
C     TEST IF R - REFINEMENT IS REQUESTED
      IF (IDOB.EQ.0) GO TO 20
      IDOB1=IDOB
      IDOB=0
      NV=NV-NVD
      IF (NCYCLE.NE.0) GO TO 20
   19 NCYCLE=NCYC+IDOB1
      NCB=NCYC
      NNCY=NCYCLE+1
      IDOB=1
      NV=NV+NVD
      IDOB1=0
      ICVTST=0
C*** REFINEMENT LOOP
   20 DO 103 NCYR=NCB,NNCY
      IF (NCYC.EQ.NCYCLE.AND.IDOB1.GT.0) GO TO 19
C*** COUNT REFINEMENT CYCLES
      NCYC=NCYC+1
      NCC=NCYC
      IF (IGIT.NE.1) GO TO 11
      CALL SPEICH
   11 CALL DISDER
      IF (IRNG.GT.0.AND.RWERT.GT.RWM) GO TO 16
      IF (IMAT.EQ.0.OR.NCYC.GT.NCYCLE) GO TO 3
      IF (IMAT.EQ.2) GO TO 2
      IF (NCYC.NE.NCYCLE) GO TO 3
    2 MATCOR=0
      CALL MATOUT
    3 IF (IDIS.EQ.1) GO TO 4
      IF (NCYC.NE.1.AND.NCYC.NE.NCYCLE+1) GO TO 5
    4 CALL DISOUT
    5 WRITE (NTOUT,1000) RWERT,NCC
      IZEIL=IZEIL+3
      WRITE (NTOUT,1001) RSIGMA
      IF (NRE.EQ.NRH) GO TO 77
      CALL RHO
      WRITE (NTOUT,2000) ROV
      IZEIL=IZEIL+1
 2000 FORMAT (43H RHO=SUM(W*(DO-D))**2 + SUM(W*COND)**2    =,F13.8)
 1000 FORMAT (//' R=SQRT( SUM(W*(DO-D))**2 / SUM(W*DO)**2 )= ',F13.8,
	1 5X,' BEFORE CYCLE',I4)
 1001 FORMAT (43H SIGMA=SQRT( SUM(W*(DO-D))**2 / (M-NV) )  =,F13.8)
C     TEST IF JOB IS FINISHED
C***  WRITE PARAMETERS ON NFILEA IF CONVERGENCE AND IF DESIRED
   77 CALL KOPF1
    7 IF (ICVTST.EQ.0.AND.NCYC.LE.NCYCLE) GO TO 16
      IF (IPCH.EQ.1) GO TO 31
      IF (ICVTST.GT.0.OR.NCYC.GT.NCYCLE) GO TO 105
      GO TO 24
   16 CALL INVERT
      CALL NEWPAR
      IF (IORT.EQ.1) GO TO 6
      IF (NCYC.NE.1.AND.NCYC.NE.NCYCLE) GO TO 28
      IF (NCYC.EQ.1.AND.IRNG.GT.0) GO TO 28
    6 CALL PAROUT
   28 IF (IPCH.NE.2) GO TO 24
   31 WRITE (NFILEA,1010) (ITIT(I),I=1,18)
 1010 FORMAT (  6HTITLE ,18A4)
      NCYP=NCYC
      IF (NCYC.GT.NCYCLE) NCYP=NCYCLE
      WRITE (NFILEA,1011) NCYP,NCAPID,KRNG
 1011 FORMAT (  6HCYCLE ,I4,2I5)
      IF (IGIT.NE.1) GO TO 27
C***  WRITE PARAMETERS ON NFILEA
      DO 32 I=1,3
   32 A1(I) = A(I+3)/PI180
      WRITE(NFILEA,1007) LSYS,(A(I),I=1,3), (A1(I),I=1,3)
 1007 FORMAT (6HCELL  ,1X,A4,9X,3F10.5,3F10.6)
   27 DO 100 I=1,N1
  100 WRITE (NFILEA,1002) ID(I),(X(K,I),K=1,3)
 1002 FORMAT (6HNATOM ,1X,A6,3F8.5)
      WRITE (NFILEA,1003) NTIN
 1003 FORMAT (6HFILES ,I3)
      IF (IDOB.EQ.0.AND.NCAPID.EQ.0) GO TO 23
      DO 101 J=1,M
      K=LJ(J)
      IF (IDOB.EQ.0) GO TO 25
      IK=LDR(K)
      IF (IK.EQ.0) GO TO 25
      SP=RAT(K)*DOB(K)
      IS=ISYMB(IK)
      GO TO 26
   25 IS=ISYMB(K)
      SP=0.0
   26 CONTINUE
      K1=NN(J,1)
      K2=NN(J,2)
      SQRTW=W(J)
      IF (MW(J).EQ.2) SQRTW=SQRTW/WFAC
      SQRTW=SQRT(SQRTW)
      WRITE (NFILEA,1004) ID(K1),ID(K2),DOB(K),SQRTW,IS,SP
 1004 FORMAT (6HDISTAN,1X,2(A6,1X),2F10.5,A1,F10.5)
  101 CONTINUE
   23 WRITE (NFILEA,1003) NTIN
   24 IF (ICVTST.GT.0.OR.NCYC.GT.NCYCLE) GO TO 105
      IF (IABS(ICOR)-1) 9,17,8
   17 IF (NCYC.NE.NCYCLE.OR.INEW.NE.0) GO TO 9
    8 IF (IDOB.EQ.1) GO TO 108
      MATCOR=1
      CALL MATOUT
      GO TO 9
  108 WRITE (NTOUT,2008)
      IZEIL=IZEIL+1
      CALL KOPF1
      ICOR=0
    9 IF (NCYC.GT.NCYCLE) GO TO 105
      IF (ICAL.NE.0) GO TO 103
C***  CONVERGENCE TEST
      DO 102 I=1,NV
      IF (ABS(DELV(I)).GE.CVGTST) GO TO 103
  102 CONTINUE
      ICVTST=1
      NCC=NCYC+1
      WRITE (NTOUT,1006) NCYC
      IZEIL=IZEIL+1
      CALL KOPF1
 1006 FORMAT (//' CONVERGENCE TEST POSITIVE AFTER CYCLE',I4/)
C***  PUT OUT APPROX. EIGENVALUES OF NORMAL MATRIX
      IF (IDOB.EQ.1) GO TO 60
      MEIG=1
      CALL MATOUT
      IF (IDOB1.GT.0) GO TO 19
      IF (INEW.NE.0) GO TO 70
   60 IF (IABS(ICOR).NE.1.OR.NCYC.EQ.NCYCLE) GO TO 70
      IF (IDOB.EQ.1) GO TO 107
      MATCOR=1
      CALL MATOUT
      GO TO 70
  107 WRITE (NTOUT,2008)
      IZEIL=IZEIL+1
      CALL KOPF1
      ICOR=0
   70 IF (IORT.NE.0.OR.NCYC.EQ.1.OR.NCYC.EQ.NCYCLE) GO TO 80
      CALL PAROUT
   80 CALL RHO
      IF (IRNG.NE.0.AND.RWERT.GE.RWM) GO TO 5
      GO TO 4
  103 CONTINUE
  105 IF (INEW.EQ.0.OR.ICOR.EQ.0) GO TO 205
      IF (IDOB.EQ.1) GO TO 109
      INEW1=INEW
      INEW=0
      IF (IGIT.EQ.1) CALL SPEICH
      CALL DISDER
      CALL INVERT
      MATCOR=1
      CALL MATOUT
      INEW=INEW1
      GO TO 205
  109 WRITE (NTOUT,2008)
 2008 FORMAT (//'*** REMEMBER: THE CORRELATION MATRIX IS NOT ',
     1'CALCULATED IN CONST. RATIO REFINEMENT CYCLES'/)
      IZEIL=IZEIL+1
      CALL KOPF1
      ICOR=0
  205 IF (RWERT.GT.RWM) GO TO 1
      IF (NCAPID.EQ.NAPID) GO TO 110
      IF (NCYC.EQ.1.AND.NCAPID.NE.0) GO TO 110
      IF (RWERT.LT.RWM) GO TO 106
      WRITE (NTOUT,1030) RWM
      GO TO 1
  106 CALL APID
      ICVTST=0
      NCYC=0
      GO TO 20
  110 ICVTST=2
      CALL PAROUT
C
C     CALCULATE AND PRINT TOTAL SHIFTS
C
C     ADAPT X,Y AND/OR Z PARAMETERS OF DLS MODEL TO INITIAL PARAMETERS
C     IF REQUESTED
C
C
c     If you don't want a powd10 file you may cut from
c     here to '*** End of powd10 part ***' or you may
c     set nfileb = 0
c
C     Write input file for powd10?
C
      if (nfileb .gt. 0) then
        WRITE(6,3000)
 3000   FORMAT(//'    Do you want to write an input file for POWD10? ',
     1  '(Y/N) ',$)
        READ(5,3001) ANS
 3001   FORMAT(A)
      endif
c
      IF(      nfileb .gt. 0
     1   .and. (ANS(1:1).EQ.'Y' .OR. ANS(1:1).EQ.'y')) THEN
        WRITE(NFILEB,3500)
 3500   FORMAT(' 0 1 1 0 0 0 0 1 0 1 0 1 0 1')
        WRITE(6,3010)
 3010   FORMAT('    Title Card 1 > ',$)
        READ(5,3001) TIT1
        WRITE(NFILEB,3510) TIT1
 3510   FORMAT(A72)
        WRITE(6,3020)
 3020   FORMAT('    Title Card 2 > ',$)
        READ(5,3001) TIT2
        WRITE(NFILEB,3510) TIT2
        WRITE(6,3030)
 3030   FORMAT('    Title Card 3 > ',$)
        READ(5,3001) TIT3
        WRITE(NFILEB,3510) TIT3
        WRITE(6,3040)
 3040   FORMAT('    Title Card 4 > ',$)
        READ(5,3001) TIT4
        WRITE(NFILEB,3510) TIT4
        WRITE(6,3041)
 3041   FORMAT('    Space Group Number > ',$)
        READ(5,*)NSG
        NTCS=0
        WRITE(6,3042)
 3042   FORMAT('    Theta Compensating Slit? (Y/N) > ',$)
        READ(5,3001) ANS
        IF(ANS(1:1).EQ.'Y' .OR. ANS(1:1).EQ.'y') NTCS=1
        WRITE(6,3043)
 3043   FORMAT(' Do you want a wavelength other than Cu K-alpha ? ',
     1  '(Y/N) ',$)
        READ(5,3001) ANS
        IF(ANS(1:1).EQ.'Y' .OR. ANS(1:1).EQ.'y') THEN
         WRITE(6,3044)
 3044    FORMAT(' Desired wavelength (in Angstroms) ',$)
         READ(5,*) XLAMDA
         WRITE(NFILEB,3515) NSG,N1,NTCS,XLAMDA
 3515    FORMAT('XX',I4,I3,'  1  0  0  0  0  0  0  0',I3,'  0  0'
     1   '  0  0  0.0   0.0',F10.5)
        ELSE
         WRITE(NFILEB,3520)NSG,N1,NTCS
 3520    FORMAT('CU',I4,I3,'  1  0  0  0  0  0  0  0',I3,'  0  0')
        END IF
        WRITE(NFILEB,3530)
 3530   FORMAT(' 21 27 31 34 36')
        DO 3532 I=1,3
 3532     A1(I) = A(I+3)/PI180
        WRITE(NFILEB,3540) (A(I),I=1,3),(A1(I),I=1,3)
 3540   FORMAT(3F10.4,3(F7.3,3X))
        DENS = 1.9
        WRITE(6,3045)
 3045   FORMAT('    Maximum 2 theta > ',$)
        READ(5,*)TTMAX
        WRITE(6,3046)
 3046   FORMAT('    Space Group Multiplicity > ',$)
        READ(5,*)SCAL
        WRITE(NFILEB,3550)DENS,TTMAX,SCAL
 3550   FORMAT(' 0 0 0   0.0 100.0',2F6.2,F6.1)
        DO 3100 I=1,N1
         REWIND 1
         WRITE(1,2990) ID(I)
 2990    FORMAT(A6)
         REWIND 1
         READ(1,2991) ATID
 2991    FORMAT(A2)
C BK     IF (KTYPE(ITYPE(I)).EQ.'SI') THEN
         IF (KTYPE(ITYPE(I)).EQ.IBKSI) THEN
           SF=4.0
           BEQ=1.0
           GO TO 3555
         ENDIF
C BK     IF (KTYPE(ITYPE(I)).EQ.'AL') THEN
         IF (KTYPE(ITYPE(I)).EQ.IBKAL) THEN
           SF=3.0
           BEQ=1.0
           GO TO 3555
         ENDIF
C BK     IF (KTYPE(ITYPE(I)).EQ.'O') THEN
         IF (KTYPE(ITYPE(I)).EQ.IBKO) THEN
           SF=1.0
           BEQ=2.0
           GO TO 3555
         ENDIF
C BK     IF (KTYPE(ITYPE(I)).EQ.'NA') THEN
         IF (KTYPE(ITYPE(I)).EQ.IBKNA) THEN
           SF=2.0
           BEQ=3.0
           GO TO 3555
         ENDIF
C BK     IF (KTYPE(ITYPE(I)).EQ.'P') THEN
         IF (KTYPE(ITYPE(I)).EQ.IBKP) THEN
           SF=5.0
           BEQ=1.0
           GO TO 3555
         ENDIF
 3555    WRITE(NFILEB,3560) ID(I),SF,XMUL(I),(X(K,I),K=1,3),BEQ
 3100    CONTINUE
 3560   FORMAT(A6,2X,2F6.1,4F10.6)
        WRITE(6,3047)
 3047   FORMAT(' Do you want to input a halfwidth other than the ',
     1  'default function ? (Y/N) ',$)
        READ(5,3001) ANS
        IF(ANS(1:1).EQ.'Y' .OR. ANS(1:1).EQ.'y') THEN
         WRITE(6,3048)
 3048    FORMAT(' Halfwidth desired ',$)
         READ(5,*) WIDTH
         WRITE(6,3049)
 3049    FORMAT(' 2 Theta value for this halfwidth ',$)
         READ(5,*) TTX
         WRITE(NFILEB,3566) WIDTH,TTX
 3566    FORMAT(2F6.2,' 100.0 0 0 0 1 0 0 0    0.10    0    0',
     1   '    0    0     0 0 0 1')
        ELSE
         WRITE(NFILEB,3570)
 3570    FORMAT('  0.00  0.00 100.0 0 0 0 1 0 0 0    0.10    0    0',
     1   '    0    0     0 0 0 1')
        END IF
 3575   WRITE(6,3050)
 3050   FORMAT('    POWD10 output file name (A6) > ',$)
        READ(5,3001)P10FN
        WRITE(NFILEB,3580)P10FN
 3580   FORMAT(A6)
        WRITE(NFILEB,3590)
 3590   FORMAT(' ')
        REWIND NFILEB
        WRITE(6,3060)
 3060   FORMAT(/
     1  '    A file called POWD10.INP has been created for you'//
     2  '    Multiplicities for atoms in special positions will have'/
     3  '    to be changed before running POWD10 with this input '/
     4  '    file.'/)
      ENDIF
c
c     *** End of powd10 part ***
C
      DO 120 K=1,3
      IF (IMOV(K).EQ.0) GO TO 120
      IADAPT=1
      WRITE (NTOUT,1025) IXYZ(K)
      IZEIL=IZEIL+2
      CALL KOPF1
      SUMDX=0.0
      DO 115 II=1,N1
  115 SUMDX=SUMDX+XINI(K,II)-X(K,II)
      DX=SUMDX/FLOAT(N1)
      DO 116 I=1,N1
  116 X(K,I)=X(K,I)+DX
  120 CONTINUE
      IF (IADAPT.EQ.0) GO TO 1
      IF (NRE.GT.0) WRITE (NTOUT,1020)
      ICVTST=3
      CALL PAROUT
      GO TO 1
 1020 FORMAT (//' *** CAUTION: THE ADAPTED PARAMETERS MAY NO MORE AG'
     1,'REE WITH THE LINRES CONDITIONS: T H I N K')
 1025 FORMAT (//' TRANSLATE  REFINED DLS MODEL ALONG 'A2, 55HTO MINIMI
     1SE AVERAGE DEVIATIONS FROM INITIAL COORDINATES)
 1030 FORMAT (//' AGREEMENT FACTOR IS GREATER THAN ',F5.2, 18H , SKIP
     1 APID CYCLE)
      END
      SUBROUTINE DATIN
C      PURPOSE... THE ONLY ROUTINE THAT READS INPUT DATA.
C       NO OTHER PART OF THE PROGRAM CONTAINS A READ-STATEMENT.
      DIMENSION XX(3), ISYST(8)
      DIMENSION FIZ(5), KOOR(5), IDOM(5), CHAR(5), KKAT(2), CHAR2(6)
      DIMENSION IX(3), KARTE(14), KONT(18), KAT(3), KREF(6), NCON(9)
      REAL*8 KONT,KARTE,LATOM,K1,CHAR1,WCODE(201),ID,KKAT,IDOM,NCON,
     1        LABEL
      REAL IX,IHD,IHI,IHL,IHR,ISM,IBL1,IDEC,IFAK,IGIN,KOOR,KREF,
     1     ISYMB,LCELL,K3
      INTEGER BT
      CHARACTER*10 KSYMB
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK4/ RN,RSIGMA,RWERT,ICVTST
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK7/ IDEC(100),KA(2,150),LCELL(6)
      COMMON /BLOCK8/ SI(3,150),SIGN(3,3,150),SYMOPS(10,200),B(3,3,200),
     1BI(3,200),IND,MTR(200)
      COMMON /BLOCK9/ DELV(300),IFELD(6),IV(300),IV1(300),KVAR(3,200)
      COMMON S(50000)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
      COMMON /BLOC12/ SYMSIG(3,3,40),KATOR(400),NN11(400,2)
      COMMON /BLOC13/ ISYS,CVGTST,LSYS
      COMMON /BLOC14/ IZ1(35),WR(35),BEDING(35),RESKO(35,20),KARI(35,20)
     1,LINH(35),LIND(35)
      COMMON /BLOC15/ NAPID,IRNG,KRNG
      COMMON /BLOC16/ ITYPE(200),BT(15),IA,NZA,ICON(40,11),MD(40,18),NCA
     1PID,CAPID(15,4),IAT,IAB,IBT(40,4),DOT,WEIGHT(15,3),DELTAX(3),DELTA
     20(3),DA,DB,ANG,TOT(40,4),OTO(40,6),PI180,IMOV(3),MSYMAB(200),DOBIN
     3(400),KTYPE(20),IDK
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
      COMMON /BLOC18/ ROV,NC,NA,NS,IVF,IRNG1
      COMMON /BLOC19/ RDD,ROTOL,ROTOU,RTOTL,RTOTU
C BK      COMMON /BLOC22/ NUMB(10),IXYZ(3),LCENT(4),ISLASH,IBLANK,ICOMMA,
C BK     1IPLUS,IMINUS,ISTAR
      COMMON /BLOC22/ IXYZ(3),LCENT(4),ISLASH,IBLANK,ICOMMA,
     1IPLUS,IMINUS,ISTAR
      COMMON /BLOC30/ ISPACE,INPF,IPRINT,ICENT,NSYM,LTRAN,RL(3,4),
     $ KSYMB(3,48),RS(48,3,4),ITRAN(30),IJUNK(30)
      COMMON /BLOC31/ AXIS(3),KEY(3),DLIM(3),XI(3),XJS(3),
     1 TJUNK(3),SEP(3),ISGA(200),NAMX(6)
     2 ,RADI(20),ITET(200)
C BK
      COMMON /BLOCKBK/ NUMB(10)
      DATA IHD,IHL/'D','L'/
      DATA KONT/'TITLE','DLS-76','CELL','ATOM','SYMEQ','LINRES','NOREF',
     1          'DISTAN','FILES','BONDIS','END','FINISH','TETCON',
     2          'FACTOR','CYCLE','NATOM','RANGES','GENER'/
      DATA ISYST/'TRIC','MON1','MON2','ORT','TET','HEX','RHO','CUB'/
      DATA LATOM/'ATOMS'/
      DATA IX/'X','Y','Z'/
      DATA I40/40/
      DATA CHAR/'R','I',' ',',','='/
      DATA CHAR1/'      '/
      DATA CHAR2/'A','B','C','ALFA','BETA','GAMA'/
      IHR=CHAR(1)
      IHI=CHAR(2)
      IBL1=CHAR(3)
      LCELL(1)=CHAR2(1)
      LCELL(2)=CHAR2(2)
      LCELL(3)=CHAR2(3)
      LCELL(4)=CHAR2(4)
      LCELL(5)=CHAR2(5)
      LCELL(6)=CHAR2(6)
C*** CHECK IF RANDOM COORDINATES DESIRED
      IF (IRNG.EQ.0) GO TO 56
      IF (KRNG.LT.IRNG) GO TO 58
      KRNG=0
      GO TO 56
   58 KRNG=KRNG+1
      IRNG1=IYY
      GO TO 57
   56 CONTINUE
      NFUNC=0
      IRNG=0
      PI180=3.141592654/180.
C
C FOLLOWING FILE STREAM SETUP MOVED TO OPENFILES. JAN 2002 RJA
C     NTIN= 7
C      NTOUT=12
C      NFILEA=8
      NTIN2=9
      NTIN3=10
c
c     If you       want a powd10 file set nfileb = 11
c     If you don't want a powd10 file set nfileb = 0
c     NFILEB=11
      nfileb = 0
c
      NREAD=NTIN
C***  SCRATCH FILE NTINF TO HANDLE INPUT DATA
      OPEN (UNIT=1,STATUS='SCRATCH')
      NTINF=1
      IA=0
      MOLD=1
      MCON=0
      ISYS=0
      M1=0
      M2=0
      NZA=0
      IERR=0
      NATYP=0
      MEQD=0
      ICYC=0
C***  DIMENSION LIMITS OF SOME ARRAYS
      IDIMLI=35
      IDIMAP=15
      IDIMNZ=40
      IDKTYP=20
      IDIMM=400
      IDIMN=200
      IDIMNV=300
      IDNEQU=150
      IDIMIZ=20
      IDIMRT=40
      IDIMCT=200
      IDK=IDKTYP+5
      CVGTST=0.0001
      WRF=1.
      WFAC=1.0
      CORR=0.5
      DOLD=0.0
      RDD=0.03
      ROTOL=104.5
      ROTOU=114.5
      RTOTL=115.5
      RTOTU=175.5
      DO 90 I=1,3
   90 FUDGE(I)=1.0
      DO 99 I=1,18
   99 ITIT(I)=CHAR(3)
      IDEC(1)=IBL1
      DO 400 I=31,41
  400 IDEC(I)=IBL1
      IZZ1=0
      NDO=0
      NVD=0
      DO 131 I=1,IDIMM
  131 ISYMB(I)=IBL1
      IZZ2=0
      N=0
      N1=0
      N2=0
      NEQU=0
      NRE=0
      NRH=0
      KF=0
      M=0
      MYS=0
      RN=0.0
      IPAGE=0
      NTIT=0
C***  READ INPUT CARD AND WRITE ON NTINF
C     STORE THE CODE NAME OF THE DATA CARD IN ARRAY WCODE
      REWIND NTINF
      ICOUNT=0
    1 READ(NREAD,1000,END=611) KARTE
      IF(KARTE(1).EQ.KONT(12)) GO TO 611
      IF(KARTE(1).EQ.KONT(15)) GO TO 622
      WRITE (NTINF,1000) KARTE
 1000 FORMAT (13A6,A2)
      ICOUNT=ICOUNT+1
      WCODE(ICOUNT)=KARTE(1)
      IF(ICOUNT.GT.IDIMCT) GO TO 11
      IF(KARTE(1).EQ.KONT(11).OR.KARTE(1).EQ.KONT(9)) GO TO 11
      GO TO 1
   11 REWIND NTINF
      DO 2000 IJ=1,ICOUNT
      DO 100 I=1,18
      IF (WCODE(IJ).EQ.KONT(I)) GO TO (601,602,603,604,605,606,607,608,
     1609,613,610,611,620,621,622,623,624,627),I
  100 CONTINUE
      WRITE (NTOUT,1001)
 1001 FORMAT (1HO, 34HTHE FOLLOWING CARD IS NOT ALLOWED:)
      WRITE (NTOUT,1040) KARTE
      IERR=IERR+1
      GO TO 2000
  140 BACKSPACE NTINF
      READ (NTINF,1000) KARTE
      WRITE (NTOUT,1040) KARTE
 1040 FORMAT ( 13H CARD IMAGE: ,13A6,A2)
      IERR=IERR+1
      GO TO 2000
C***  TITLE CARD
  601 READ (NTINF,1002) ITIT
 1002 FORMAT (6X,18A4)
      IF (NTIT.GT.0) GO TO 2000
      CALL KOPF
      NTIT=1
      GO TO 2000
C***  DLS-76 CARD
  602 IZZ1=1
      READ (NTINF,1003) INEW,NDIAG,IRNG,NCYCLE,IDOB,NAPID,IGIT,NC,NA,NS,
     1 ICAL,IVF,IMOV(1),IMOV(2),IMOV(3),IORT,IDIS,IMAT,ICOR,IPCH
 1003 FORMAT (7X,I3,14I2,2X,5I2)
      NRAGA=0
      IF (INEW.EQ.0.OR.IABS(ICOR).NE.2) GO TO 670
      WRITE (NTOUT,1060)
 1060 FORMAT (' *** REMEMBER: WHEN USING NEWTON-RAPHSON REFINEMENT',
     1' THE CORRELATION MATRIX IS ONLY CALCULATED AFTER THE FINAL',
     2' CYCLE')
      IZEIL=IZEIL+2
      ICOR=ICOR/2
  670 IF (INEW.LT.0) NRAGA=-INEW
      IF (NAPID.EQ.0) GO TO 2000
      IF (IDOB.EQ.0) GO TO 2000
      IDOB=0
      WRITE (NTOUT,2010)
 2010 FORMAT (//'***   REFINEMENT OF PRESCRIBED DISTANCES IS NOT PO',
     1'SSIBLE IN AN APID-RUN')
      GO TO 2000
C***  CELL CARD
  603 IZZ2=1
      READ (NTINF,1004) LSYS,A
 1004 FORMAT (7X,A4,9X,6F10.4)
      DO 97 K=4,6
      IF (A(K).LT.1.0) GO TO 96
      A(K)=A(K)*PI180
      GO TO 97
   96 A(K)=ACOS(A(K))
   97 CONTINUE
      DO 101 I=1,6
      IF (IGIT.NE.0) GO TO 61
      IGIN(I)=IHI
      GO TO 101
   61 IGIN(I)=IHR
  101 CONTINUE
      ISYS=1
      DO 98 K=1,8
      IF (LSYS.EQ.ISYST(K)) GO TO (70,29,30,31,32,32,33,34), K
   98 CONTINUE
      WRITE (NTOUT,2045) LSYS
 2045 FORMAT (//' ERROR ON CELL CARD: ILLEGAL CRYSTAL SYSTEM   ',
     1' INDICATOR ',A4)
      IERR=IERR+1
      GO TO 70
C
C***  MONOCLINIC FIRST SETTING
C
   29 IGIN(4)=IHI
      IGIN(5)=IHI
      GO TO 70
C
C***  MONOCLINIC SECOND SETTING
C
   30 IGIN(4)=IHI
      IGIN(6)=IHI
      GO TO 70
C
C***  ORTHORHOMBIC
C
   31 IGIN(4)=IHI
      IGIN(5)=IHI
      IGIN(6)=IHI
      GO TO 70
C***   TRIGONAL,HEXAGONAL,TETRAGONAL
   32 IGIN(2)=IHD
      IGIN(4)=IHI
      IGIN(5)=IHI
      IGIN(6)=IHI
      A(2)=A(1)
      IF (LSYS.EQ.ISYST(6)) A(6)=120.*PI180
      ISYS=2
      GO TO 70
C***   RHOMBOHEDRAL
   33 IGIN(2)=IHD
      IGIN(3)=IHD
      IGIN(5)=IHD
      IGIN(6)=IHD
      A(2)=A(1)
      A(3)=A(1)
      A(5)=A(4)
      A(6)=A(4)
      ISYS=3
      GO TO 70
C***   CUBIC
   34 IGIN(2)=IHD
      IGIN(3)=IHD
      IGIN(4)=IHI
      IGIN(5)=IHI
      IGIN(6)=IHI
      A(2)=A(1)
      A(3)=A(1)
      ISYS=4
   70 CALL SPEICH
      GO TO 2000
C***  ATOM CARD
  604 N1=N1+1
      N=N1
      IF (N.LE.IDIMN) GO TO 75
      WRITE (NTOUT,2090) IDIMN
 2090 FORMAT (' DIMENSION CHECK: NUMBER OF ATOM AND SYMEQ', 
	1' CARDS LIMITED TO',I4)
      IERR=IERR+1
      N1=N1-1
      N=N-1
      GO TO 2000
   75 MSYMAB(N1)=N1
      DO 102 K=1,3
  102 IFAK(K,N1)=IHR
      SYMOPS(9,N1)=IBL1
      SYMOPS(10,N1)=IBL1
      READ (NTINF,1005) ID(N1),(X(K,N1),K=1,3),ITYPE(N1),
     1(IDEC(I),I=2,31),XMUL(N1)
 1005 FORMAT (7X,A6,3F8.5,2X,A3,1X,30A1,F6.1)
      IF (XMUL(N1).LT.0.0005) XMUL(N1) = 1.0
      IF (NATYP.NE.0) GO TO 300
      NAT=1
      NATYP=1
      GO TO 310
  300 DO 305 NT=1,NATYP
      NAT=NT
      IF (KTYPE(NT).EQ.ITYPE(N1)) GO TO 315
  305 CONTINUE
      NAT=NAT+1
      NATYP=NAT
      IF (NATYP.LE.IDKTYP) GO TO 310
      WRITE (NTOUT,2005) IDKTYP
 2005 FORMAT (//'  DIMENSION CHECK: NUMBER OF ATOM TYPES LIMITED',
     1' TO ',I3)
      NATYP=IDKTYP
      IERR=IERR+1
C***  INCREASE DIMENSION OF KTYPE AND SET IDKTYP TO NEW VALUE
      GO TO 316
  310 KTYPE(NAT)=ITYPE(N1)
  315 ITYPE(N1)=NAT
  316 DO 401 I=2,31
      IF (IDEC(I).NE.IBL1) GO TO 402
  401 CONTINUE
      MTR(N1)=0
      DO 500 I=1,3
      BI(I,N1)=0.0
      DO 500 K=1,3
  500 B(I,K,N1)=0.0
      DO 501 I=1,3
  501 B(I,I,N1)=1.0
      GO TO 407
  402 MTR(N1)=N1
      DO 403 K=1,3
      DO 404 I=2,31
      IF (IDEC(I).NE.IX(K)) GO TO 404
      GO TO 403
  404 CONTINUE
      IFAK(K,N1)=IHI
  403 CONTINUE
      K=1
      DO 502 I=2,31
      IF (IDEC(I).EQ.CHAR(4)) GO TO 503
      IF (IFAK(K,N1).NE.IHR) GO TO 502
      IF (IDEC(I).NE.IBL1.AND.IDEC(I).NE.IX(K)) GO TO 504
      GO TO 502
  503 K=K+1
  502 CONTINUE
      GO TO 505
  504 WRITE (NTOUT,2085) ID(N1),(IDEC(I),I=2,31)
 2085 FORMAT (' ERROR IN SPECIAL POSITION OF ',A6,
	1' INDEPENDENT, COOR. INCORRECTLY INTRODUCED: ',30A1)
      IERR=IERR+1
      GO TO 2000
  505 BACKSPACE NTINF
      READ (NTINF,1105) (SYMOPS(I,N1),I=1,8)
 1105 FORMAT (43X,7A4,A2)
      IND=1
      CALL SYMOP
      DO 405 I=1,3
      IF (IFAK(I,N1).EQ.IHR) GO TO 405
      DO 406 K=1,3
      IF (ABS(B(I,K,N1)).LT.0.1) GO TO 406
      IFAK(I,N1)=IHD
      GO TO 405
  406 CONTINUE
  405 CONTINUE
  407 DO 408 K=1,3
  408 XINI(K,N1)=X(K,N1)
      GO TO 2000
C
C***  NATOM CARD  (NEW ATOM PARAMETERS)
C     CONTAINS REFINED PARAMETERS FROM PREVIOUS JOB
C
  623 READ (NTINF,1005) LABEL
      DO 410 I=1,N1
      IF (ID(I).EQ.LABEL) GO TO 415
  410 CONTINUE
      WRITE (NTOUT,2075) LABEL
 2075 FORMAT (//'*** ERROR: ATOM LABEL ',A6, 20H OF PREVIOUS JOB WAS
     1, 44H NOT ENCOUNTERED ON AN ATOM CARD IN THIS JOB)
      IERR=IERR+1
      GO TO 2000
  415 BACKSPACE NTINF
      READ (NTINF,1005) LABEL,(X(K,I),K=1,3)
      GO TO 2000
C***  SYMEQ CARD
  605 NEQU=NEQU+1
      IF (NEQU.LE.IDNEQU) GO TO 420
      WRITE (NTOUT,3000) IDNEQU
 3000 FORMAT (//'*** DIMENSION CHECK: NUMBER OF SYMEQ CARDS',
     1'LIMITED TO ',I3)
      IERR=IERR+1
      NEQU=NEQU-1
      GO TO 2000
  420 READ (NTINF,1006) KKAT,(IDEC(I),I=2,41)
 1006 FORMAT (7X,A6,1X,A6,1X,40A1)
      N2=N2+1
      N=N1+N2
      IF (N.LE.IDIMN) GO TO 80
      WRITE (NTOUT,2090) IDIMN
      IERR=IERR+1
      N2=N2-1
      N=N-1
      GO TO 2000
   80 BACKSPACE NTINF
      READ (NTINF,1030) (SYMOPS(I,N),I=1,10)
 1030 FORMAT (21X,10A4)
      ID(N)=KKAT(2)
      KA(2,NEQU)=N
      DO 103 I=1,N1
      IF (KKAT(1).NE.ID(I)) GO TO 103
      KA(1,NEQU)=I
      GO TO 2
  103 CONTINUE
      WRITE (NTOUT,1007) WCODE(IJ),NEQU,KKAT(1)
 1007 FORMAT (//' *** ERROR ON ',A6,  8H CARD NO,I3, 21H: UNKNOWN ATO
     1M LABEL ,A6)
      GO TO 140
    2 MSYMAB(N)=KA(1,NEQU)
      MTR(N)=N
      IND=2
      CALL SYMOP
      GO TO 2000
C***   LINRES CARD
  606 IF (KF.NE.0) GO TO 320
      NRE=NRE+1
      IF (NRE.LE.IDIMLI) GO TO 125
      WRITE (NTOUT,2070) IDIMLI
 2070 FORMAT (//' *** DIMENSION CHECK: NUMBER OF LINRES CARDS ',
	1	'LIMITED TO',I3)
      NRE=NRE-1
      IERR=IERR+1
      GO TO 2000
  125 IZ1(NRE)=0
      LINH(NRE)=0
      LIND(NRE)=0
      READ (NTINF,1025) WR(NRE)
 1025 FORMAT (7X,F8.1)
      IF (WR(NRE).GT.0.0) GO TO 46
C     IT IS A HARD CONSTRAINT
      LINH(NRE)=1
      NRH=NRH+1
      GO TO 320
   46 WR(NRE)=WR(NRE)*WRF
      WR(NRE)=WR(NRE)*WR(NRE)
  320 BACKSPACE NTINF
      READ (NTINF,1026) (FIZ(I),KOOR(I),IDOM(I),I=1,5),K3
 1026 FORMAT (16X,5(F4.0,1X,A1,A6),1X,A1)
      DO 123 I=1,5
      IF ((IDOM(I).EQ.CHAR1).AND.(KOOR(I).EQ.IBL1)) GO TO 123
      IZ1(NRE)=IZ1(NRE)+1
      NH=IZ1(NRE)
      IF (NH.LE.IDIMIZ) GO TO 323
      WRITE (NTOUT,3005) NRE,IDIMIZ
 3005 FORMAT (//'*** DIMENSION CHECK: NUMBER OF TERMS IN LINEAR', 
	1' RESTRICTION NR',I3, ' EXCEEDS LIMIT OF',I3)
      IERR=IERR+1
      IF (K3.EQ.CHAR(5)) GO TO 322
      KF=1
      GO TO 2000
  322 KF=0
      NRE=NRE-1
      GO TO 2000
  323 IF (FIZ(I).EQ.0.0) FIZ(I)=1.0
      RESKO(NRE,NH)=FIZ(I)
      DO 121 L=1,3
      IF (KOOR(I).NE.IX(L)) GO TO 121
      L1=L
      GO TO 49
  121 CONTINUE
      WRITE (NTOUT,1027) NRE
 1027 FORMAT (4HO***, 24H ERROR ON LINRES CARD NO,I3, 18H: MISSING X,Y,O
     1R Z)
      GO TO 140
   49 DO 122 L=1,N1
      IF (IDOM(I).NE.ID(L)) GO TO 122
      L2=L
      IF (IFAK(L1,L2).EQ.IHR) GO TO 50
      WRITE (NTOUT,1015) NRE,IX(L1),ID(L2)
 1015 FORMAT (//'*** ERROR ON LINRES CARD NO',I3, 14H : COORDINATE ,A
     11,  9H OF ATOM ,A6, 34H IS ALREADY DEPENDENT OR INVARIANT)
      GO TO 140
  122 CONTINUE
      WRITE (NTOUT,1007) WCODE(IJ),NRE,IDOM(I)
      GO TO 140
   50 KARI(NRE,NH)=3*(L2-1)+L1
  123 CONTINUE
      IF (K3.EQ.CHAR(5)) GO TO 52
      KF=1
      GO TO 2000
   52 KF=0
      IF (LINH(NRE).EQ.0) GO TO 2000
      IF (IZ1(NRE).GT.1) GO TO 53
      WRITE (NTOUT,1017) NRE
 1017 FORMAT (//'*** ERROR ON LINRES CARD NO ',I3, 
	1'FOR A HARD CONSTRAINT ON A SINGLE COORDINATE THE NOREF CARD ',
     2' MUST BE USED')
      GO TO 140
   53 GO TO 2000
C***  NOREF CARD
  607 READ (NTINF,1008) K1
 1008 FORMAT (7X,A6)
      IF (K1.EQ.KONT(3)) GO TO 701
      GO TO 702
C***  NOREF CELL
  701 IF (IGIT.NE.0) GO TO 3
      GO TO 2000
    3 BACKSPACE NTINF
      READ (NTINF,1010) KREF
 1010 FORMAT (13X,6A4)
      DO 105 I=1,6
      DO 105 K=1,6
      IF (KREF(I).NE.LCELL(K)) GO TO 105
      IGIN(K)=IHI
  105 CONTINUE
      GO TO 2000
C***  NOREF ATOM
  702 BACKSPACE NTINF
      READ (NTINF,1011) K1,(KREF(I),I=1,3)
 1011 FORMAT (7X,A6,3(A1,2X))
      IF (K1.NE.LATOM) GO TO 35
C***  ALL ATOMS INVARIANT
      DO 120 I=1,N1
      DO 120 K=1,3
  120 IFAK(K,I)=IHI
      GO TO 2000
   35 DO 106 I=1,N1
      IF (K1.NE.ID(I)) GO TO 106
      NH=I
      GO TO 4
  106 CONTINUE
      WRITE (NTOUT,1012) K1
 1012 FORMAT (//'***', 41H ERROR ON NOREF CARD: UNKNOWN ATOM LABEL ,A6)
      GO TO 140
    4 DO 107 I=1,3
      IF (KREF(I).NE.IBL1) GO TO 5
  107 CONTINUE
      DO 108 I=1,3
  108 IFAK(I,NH)=IHI
      GO TO 2000
    5 CONTINUE
      DO 109 I=1,3
      DO 109 K=1,3
      IF (KREF(I).NE.IX(K)) GO TO 109
      IF (IFAK(K,NH).EQ.IHR) GO TO 150
      WRITE (NTOUT,1014) IX(K),ID(NH)
 1014 FORMAT (//'***ERROR ON NOREF CARD: ',A1,' OF ATOM ',A6,
     1' IS ALREADY SET DEPENDENT OR INVARIANT')
      GO TO 140
  150 IFAK(K,NH)=IHI
  109 CONTINUE
      GO TO 2000
C
C     BONDIS CARD  (BOND-TYPE AND DISTANCE FUNCTION)
C
  613 IA=IA+1
      NFUNC=1
      IF (IA.LE.IDIMAP) GO TO 614
      WRITE (NTOUT,2011) IDIMAP
 2011 FORMAT (//'***', 51H DIMENSION CHECK: NUMBER OF BONDIS CARDS LIMITE
     1D TO,I5)
      IERR=IERR+1
      IA=IA-1
      GO TO 2000
  614 READ (NTINF,2015) KAT,(CAPID(IA,K),K=1,4),(WEIGHT(IA,K),K=1,3)
 2015 FORMAT (7X,3(A3,1X),1X,4E10.4,3F5.0)
C     CALCULATE BOND TYPE INDEX
      DO 330 I=1,3
      DO 321 NT=1,NATYP
      NAT=NT
      IF (KAT(I).EQ.KTYPE(NT)) GO TO 325
  321 CONTINUE
      WRITE (NTOUT,2016) KAT(I),IA
 2016 FORMAT (//'***', 20H WARNING: ATOM TYPE ,A3, 26H PUNCHED ON BONDIS
     1CARD NO,I3, 32H DOES NOT APPEAR ON AN ATOM CARD)
      IZEIL=IZEIL+2
      KAT(I)=0
      GO TO 330
  325 KAT(I)=NAT
  330 CONTINUE
      BT(IA)=KAT(1)+(KAT(2)+2)*IDK+(KAT(3)+2)*IDK*IDK
      IF (CAPID(IA,4).EQ.0.) CAPID(IA,4)=145.
      GO TO 2000
C
C***  RANGES CARD
C
  624 READ (NTINF,2017) RDD,ROTOL,ROTOU,RTOTL,RTOTU
 2017 FORMAT (10X,5F5.2)
      E=1.E-10
      IF (RDD.LE.E) RDD=0.03
      IF (ROTOL.LE.E) ROTOL=104.5
      IF (ROTOU.LE.E) ROTOU=114.5
      IF (RTOTL.LE.E) RTOTL=115.5
      IF (RTOTU.LE.E) RTOTU=175.5
      GO TO 2000
C
C***  TETCON CARD
C
  620 MCON=1
      NZA=NZA+1
      IF (NZA.LE.IDIMNZ) GO TO 805
      WRITE (NTOUT,2025) IDIMNZ
 2025 FORMAT (//'***', 51H DIMENSION CHECK: NUMBER OF TETCON CARDS LIMITE
     1D TO,I4)
      IERR=IERR+1
      NZA=NZA-1
      GO TO 2000
  805 READ (NTINF,1020) NCON
 1020 FORMAT (13X,9(A6,1X))
C     IDENTIFY LABELS
C
C     TEST IF CENTRAL ATOM IS ON AN ATOM OR SYMEQ CARD
      DO 806 I=2,11
  806 ICON(NZA,I)=0
      DO 810 I=1,N
      IF (NCON(1).EQ.ID(I)) GO TO 815
  810 CONTINUE
      WRITE (NTOUT,2030) NZA
 2030 FORMAT (//'***ERROR ON TETCON CARD NO', I3, 66H: CENTRAL ATOM
     1 NOT SUPPLIED ON AN ATOM OR SYMEQ CARD, CARD IGNORED)
      NZA=NZA-1
      GO TO 140
  815 ICON(NZA,1)=I
      NZAT=NZA-1
      IF (NZAT.LE.0) GO TO 816
      DO 811 NZ=1,NZAT
      IF (ICON(NZA,1).NE.ICON(NZ,1)) GO TO 811
      WRITE (NTOUT,2060) NCON(1),NZ,NZA
 2060 FORMAT (//'***', 13H ERROR: ATOM , A6, 30H APPEARS TWICE AS CENTRAL
     1 ATOM, 33H ON A TETCON CARD (TETCON CARD NO,I3,  4H AND,I3,  1H)/5
     2X, 12HCARD IGNORED)
      NZA=NZA-1
      GO TO 140
  811 CONTINUE
C     CHECK REMAINING LABELS
  816 DO 830 II=2,9
      IF (NCON(II).EQ.CHAR1) GO TO 830
      DO 820 I=1,N
      IF (NCON(II).EQ.ID(I)) GO TO 825
  820 CONTINUE
      WRITE (NTOUT,2035) NZA,NCON(II),(NCON(L),L=1,9)
 2035 FORMAT (//'***', 24H ERROR ON TETCON CARD NO,I3, 16H: UNKNOWN LABEL
     1 ,A6/5X, 12HCARD IMAGE: ,9(A6,1X))
      IERR=IERR+1
      GO TO 830
  825 ICON(NZA,II)=I
  830 CONTINUE
      GO TO 2000
C***  DISTANCES
  608 IF (MCON.NE.1) GO TO 665
      IF (NFUNC.GT.0) GO TO 655
      WRITE (NTOUT,2055)
      IERR=IERR+1
      GO TO 660
  655 CALL SETUP
      M1=M
      NDO = M
  660 MCON=2
  665 M2=M2+1
      M=M1+M2
      IF (M.LE.IDIMM) GO TO 666
      WRITE (NTOUT,2036) IDIMM
 2036 FORMAT (//'***', 40H DIMENSION CHECK: NUMBER OF DISTAN CARDS, 11H L
     1IMITED TO,I4)
      IERR=IERR+1
      M2=M2-1
      GO TO 2000
  666 READ (NTINF,1016) KKAT,D1,W1,ISM,RTT
      IF (W1.LE.0.0.AND.M2.EQ.1) W1=1.0
 1016 FORMAT (7X,2(A6,1X),2F10.5,A1,F10.5)
      DO 905 L=1,2
      DO 900 I=1,N
      IF (KKAT(L).NE.ID(I)) GO TO 900
      NN(M,L)=I
      GO TO 905
  900 CONTINUE
      WRITE (NTOUT,1007) WCODE(IJ),M2,KKAT(L)
      M2=M2-1
      GO TO 140
  905 CONTINUE
C
C     TEST IF NEW DISTANCE IS EQUIVALENT TO PREVIOUS ONE
C
      NN11(M,1)=MSYMAB(NN(M,1))
      NN11(M,2)=MSYMAB(NN(M,2))
      IF (M.EQ.1) GO TO 930
      MM1=M-1
      DO 925 J=1,MM1
      JS=J
      MS=M
      CALL EQUI
      IF (ITEST.EQ.0) GO TO 925
      IF (J.GT.M1) GO TO 920
      IF (ICYC.EQ.1) GO TO 907
      MEQD=1
  907 IF (D1.GT.0.0) GO TO 910
      IF (DOLD.GT.0.0) GO TO 908
      WRITE (NTOUT,1013)
      IERR=IERR+1
      MOLD=JS
      D1=DOB(JS)
      W1=W(JS)
      GO TO 915
  908 D1=DOLD
  910 DOB(JS)=D1
  915 DOBIN(JS)=D1
      DOLD=D1
      W(JS)=W(MOLD)
      IF (W1.GT.0.0) W(JS)=W1
      MOLD=JS
      M2=M2-1
      GO TO 2000
  920 IF (ICYC.EQ.1) GO TO 922
      WRITE (NTOUT,2001) ID(NN(M,1)),ID(NN(M,2)),ID(NN(JS,1)),ID(NN(JS,2
     1))
 2001 FORMAT (//' DISTANCE ',A6,  2H- ,A6, 31H IS SYMMETRICALLY EQUIVALE
     1NT TO, 10H DISTANCE ,A6,  2H- ,A6, 34H ALREADY SPECIFIED. IT IS IG
     2NORED.)
      IZEIL=IZEIL+1
      CALL KOPF1
  922 M2=M2-1
      MOLD=JS
      IF (D1.GT.0.0) DOLD=D1
      GO TO 2000
  925 CONTINUE
  930 IF (D1.GT.0.0) GO TO 37
      IF (M2.GT.1) GO TO 36
      WRITE (NTOUT,1013)
 1013 FORMAT (//'***', 41H ERROR ON DISTAN CARD NO 1: MISSING FIRST,
     120H PRESCRIBED DISTANCE)
      IERR=IERR+1
      GO TO 37
   36 IF (W1.LE.0.0) W(M)=W(MOLD)
      LJ(M)=LJ(MOLD)
      GO TO 43
   37 IF (W1.LE.0.0) GO TO 54
      W(M)=W1
      GO TO 55
   54 W(M)=1.0
   55 IF (IDOB.EQ.0) GO TO 38
      IF (NDO.EQ.0) GO TO 38
      DO 110 K=1,NDO
      HILF=ABS(D1-DOB(K))
      IF (HILF.GT.1.0E-06) GO TO 110
      LJ(M)=K
      GO TO 43
  110 CONTINUE
   38 NDO=NDO+1
      DOB(NDO)=D1
      DOBIN(NDO)=D1
      DOLD=D1
      LJ(M)=NDO
      IF (IDOB.EQ.0) GO TO 45
      IF (NDO.LE.IDIMRT) GO TO 60
      WRITE (NTOUT,3010) IDIMRT
 3010 FORMAT (//'***', 38H DIMENSION CHECK: NUMBER OF PRESCRIBED, 18H RAT
     1IOS LIMITED TO,I3)
      IERR=IERR+1
      NDO=NDO-1
      GO TO 2000
   60 IF (ISM.NE.IBL1) GO TO 39
      LDR(NDO)=0
      GO TO 43
   39 IF (RTT.GT.0.0) GO TO 40
      RAT(NDO)=1.0
      GO TO 41
   40 RAT(NDO)=D1/RTT
   41 IF (NVD.EQ.0) GO TO 42
      DO 111 K=1,NVD
      IF (ISM.NE.ISYMB(K)) GO TO 111
      LDR(NDO)=K
      GO TO 43
  111 CONTINUE
   42 NVD=NVD+1
      ISYMB(NVD)=ISM
      LDR(NDO)=NVD
      GO TO 43
   45 ISYMB(NDO)=ISM
   43 RN=RN+W(M)*W(M)*DOB(LJ(M))*DOB(LJ(M))
      MEFF=MEFF+1
      MW(M)=1
      MOLD=M
      GO TO 2000
C
C     GENER CARD
C
  627 MCON = 1
      READ(NTINF,2100) ISPACE,INPF,IPRINT
C     DECODE(19,2100,KARTE(2))  ISPACE,INPF,IPRINT
 2100 FORMAT (10X,3I5)
      INPF = 10
C     IF(INPF.LE.0)  INPF = 8
      INSAVE = NTIN
      NTIN = INPF
      CALL GENER
      NTIN = INSAVE
      GO TO 2000
C
C*** FACTOR CARD
C
  621 READ (NTINF,2050) FUDGE,IRNG1,CVGTST,WRF,WFAC,CORR
 2050 FORMAT (10X,3F5.2,I5,4E10.4)
      IF (CVGTST.LE.1.E-14) CVGTST=0.0001
      IF (WRF.LE.1.E-10) WRF=1.0
      WFAC=WFAC*WFAC
      IF (WFAC.LE.1.E-10) WFAC=1.0
      IF (CORR.LE.1.E-10) CORR=0.5
      DD=1.E-14
      DO 91 I=1,3
      IF (FUDGE(I).LE.DD) FUDGE(I)=1.0
   91 CONTINUE
      GO TO 2000
C***  FILE ASSIGNMENT
  609 READ (NTINF,1018) NREAD,NFILES
 1018 FORMAT (6X,I3,3X,I3)
      IF (NREAD.EQ.0) NREAD=NTIN
      IF (((NREAD.EQ.NTIN).OR.(NREAD.EQ.NTIN2)).OR.(NREAD.EQ.NTIN3)) GO
     1TO 8
      REWIND NREAD
    8 IF (NFILES.EQ.0) GO TO 2000
      NFILEA=NFILES
 2000 CONTINUE
      REWIND NTINF
      ICOUNT=0
      GO TO 1
C
C***  CYCLE CARD
C
  622 ICYC=1
  729 BACKSPACE NREAD
      READ (NREAD,1045) LC,LA,LS
 1045 FORMAT (6X,I4,2I5)
      IF (LC.EQ.NC.AND.LA.EQ.NA.AND.LS.EQ.NS) GO TO 1
C***  SKIP PARAMETERS ON NREAD TILL DESIRED CYCLE IS FOUND
  730 READ (NREAD,1000,END=740) (KARTE(I),I=1,3)
  735 IF (KARTE(1).EQ.KONT(15)) GO TO 729
      GO TO 730
  740 WRITE(NTOUT,1050) NC,NA,NS,NREAD
 1050 FORMAT (//'***', 27H ERROR: PARAMETERS OF CYCLE,I3, 15H OF APID RUN
     1 NO,I3,  8H AND SET,I3, 16H ARE NOT ON FILE,I3)
      IERR=IERR+1
      NREAD=NTIN
      GO TO 1
C***  END OF DATA INPUT
  610 CONTINUE
      IF (NRE.EQ.0) GO TO 711
C     CHECK IF THE SAME VARIABLE APPEARS IN MORE THAN
C     ONE LINRES EQUATION
      KWRT=0
      NREE=NRE-1
      DO 741 L=1,NREE
      LU=L+1
      KZ1=IZ1(L)
      DO 741 I=1,KZ1
      DO 741 LL=LU,NRE
      KZ2=IZ1(LL)
      DO 742 II=1,KZ2
      IF (IABS(KARI(L,I)).NE.IABS(KARI(LL,II))) GO TO 742
      L1=MOD(IABS(KARI(L,I)),3)
      IF (L1.EQ.0) L1=3
      L2=1+(IABS(KARI(L,I))-L1)/3
      KARI(L,I)=-IABS(KARI(L,I))
      KARI(LL,II)=-IABS(KARI(LL,II))
      WRITE (NTOUT,2023) IX(L1),ID(L2),L,LL
 2023 FORMAT (//'***', 12H ATTENTION: ,A1,  9H OF ATOM ,A6, 11H APPEARS I
     1N, 19H LINRES EQUATION NO,I3,  4H AND,I3)
      IZEIL=IZEIL+1
      KWRT=1
  742 CONTINUE
  741 CONTINUE
      IF (KWRT.EQ.0) GO TO 743
      WRITE (NTOUT,2024)
 2024 FORMAT (//'   (THIS IS USUALLY ONLY POSSIBLE FOR PSEUDOHEX
     1AGONAL SYMMETRIES)')
      CALL KOPF1
  743 IF (NRH.EQ.0) GO TO 744
C     FIND VARIABLES SUITABLE FOR ELIMINATING
      DO 745 L=1,NRE
      IF (LINH(L).EQ.0) GO TO 745
      KZ1=IZ1(L)
      DO 746 I=1,KZ1
      IF (KARI(L,I).LT.0) GO TO 746
      L1=MOD(KARI(L,I),3)
      IF (L1.EQ.0) L1=3
      L2=1+(KARI(L,I)-L1)/3
      IFAK(L1,L2)=IHL
      LIND(L)=I
      GO TO 745
  746 CONTINUE
      WRITE (NTOUT,2026) L
 2026 FORMAT (//'***', 37H ERROR: CANNOT ELIMINATE ANY VARIABLE, 26H IN L
     1INEAR RESTRICTION NO ,I3)
      IERR=IERR+1
  745 CONTINUE
  744 DO 747 L=1,NRE
      KZ1=IZ1(L)
      DO 747 I=1,KZ1
      IF (KARI(L,I).LT.0) KARI(L,I)=-KARI(L,I)
  747 CONTINUE
  711 IF (MCON.NE.1) GO TO 650
      IF (NFUNC.GT.0) GO TO 645
      WRITE (NTOUT,2055)
 2055 FORMAT (//'***', 32H ERROR: NO BONDIS CARDS SUPPLIED)
      IERR=IERR+1
      GO TO 710
  645 CALL SETUP
      M1=M
C***  GENERATE RANDOM COORDINATES IF DESIRED
  650 IF (IRNG.EQ.0) GO TO 612
      IYY=IRNG1
      IF (IYY.EQ.0) IYY=12396758
      call rmarin(mod(iyy, 31329), mod(iyy, 30082))
c
   57 DO 215 I=1,N1
      DO 115 K=1,3
C***  IF THIS COORDINATE IS INVARIANT, IT IS NOT RANDOM
      IF (IFAK(K,I).NE.IHR) GO TO 115
C*** THE FOLLOWING 3 STATEMENTS GENERATE ONE RANDOM NUMBER
C    BETWEEN 0 AND 1. THE IBM RANDOM NUMBER GENERATOR RANDU(IYY,IY,XKI)
C    IS USED.  REPLACE WITH APPROPRIATE FUNCTION ON ANOTHER COMPUTER.
C*** IN THIS CASE NOTE ...
C    IRNG1 IS READ IN ON THE FACTOR CARD AND IS USED AS FIRST STARTING
C    NUMBER. (THE IBM RANDU FUNCTION REQUIRES A STARTING VALUE IYY.)
C    IF IRNG1 IS SET BY DEFAULT TO 0 IN THE INPUT, THE STARTING
C    NUMBER FOR THE RANDU FUNCTION IS ARBITRARILY SET TO 12396758.
C    XKI = GENERATED RANDOM NUMBER
C    X(K,I)= COORDINATE TO BE REPLACED BY A RANDOM NUMBER
C
C
c     CALL RANDU(IYY,IY,XKI)
c     IYY=IY
c     For information on rmarin & ranmar see below
      call ranmar(xki, 1)
C
      X(K,I)=XKI
  115 CONTINUE
      IF (NRE.NE.0) GO TO 215
      IF (MTR(I).EQ.0) GO TO 215
C     IF ATOM ON SPECIAL POSITION
      DO 216 K=1,3
  216 XX(K)=X(K,I)
      DO 217 K=1,3
      XK=0.0
      DO 218 L=1,3
  218 XK=XK+B(K,L,I)*XX(L)
  217 X(K,I)=XK+BI(K,I)
  215 CONTINUE
      IF (NRE.NE.0) GO TO 117
C***  GENERATE THE SYMMETRICALLY EQUIVALENT TO THE RANDOM ATOMS
      DO 116 I=1,NEQU
      K4=KA(1,I)
      K2=KA(2,I)
      DO 116 K=1,3
      X(K,K2)=SIGN(K,1,I)*X(1,K4)+SIGN(K,2,I)*X(2,K4)+SIGN(K,3,I)*X(3,K4
     1)+SI(K,I)
  116 CONTINUE
      IF (KRNG.GT.0.AND.NRE.EQ.0) RETURN
  117 IF (KRNG.GT.0) GO TO 716
      KRNG=1
C
C*** END RANDOM ATOMS
C
  612 CONTINUE
      M=M1+M2
      DO 210 MM=1,M
C     SQUARE WEIGHTS
      W(MM)=W(MM)*W(MM)
      L1=NN(MM,1)
      L2=NN(MM,2)
C***  STORE INDICES BEING USED IN SUBROUTINE  *DISDER*
      K4=NN11(MM,1)
      K2=NN11(MM,2)
      IF(MTR(L1).EQ.0.AND.MTR(L2).EQ.0)
     1GO TO 20
      IF (K4.EQ.K2) GO TO 21
      KATOR(MM)=-1
      GO TO 210
   20 KATOR(MM)=0
      GO TO 210
   21 MYS=MYS+1
      IF (MYS.LE.I40) GO TO 25
      WRITE (NTOUT,1022)
 1022 FORMAT (//'***', 50HDIMENSION CHECK: INCREASE DIMENSION OF SYMSIG (
     1SEE, 18H SUBROUTINE DATIN)/)
      IERR=IERR+1
      GO TO 710
   25 KATOR(MM)=MYS
      L3=MTR(L1)
      L4=MTR(L2)
      DO 200 I=1,3
      DO 200 K=1,3
      IF (L3.NE.0) GO TO 26
      SL3=FLOAT(KRONEK(I,K))
      GO TO 28
   26 SL3=B(I,K,L3)
   28 IF (L4.NE.0) GO TO 27
      SL4=FLOAT(KRONEK(I,K))
      GO TO 127
   27 SL4=B(I,K,L4)
  127 SYMSIG(I,K,MYS)=SL3-SL4
  200 CONTINUE
  210 CONTINUE
  710 CALL DEFVAR
      IF (NV.LE.IDIMNV) GO TO 709
      WRITE (NTOUT,2095) IDIMNV
 2095 FORMAT (//'***', 37H DIMENSION CHECK: NUMBER OF VARIABLES, 11H LIMI
     1TED TO,I4)
      IERR=IERR+1
  709 IF (M2.GT.0.OR.IDOB.EQ.0) GO TO 712
      WRITE (NTOUT,2039)
 2039 FORMAT (//' NOTE: THERE WERE NO DISTANCE CARDS SUPPLIED. ',
     1'CONSTANT RATIO REFINEMENT IS NOT POSSIBLE')
      IDOB=0
      IZEIL=IZEIL+2
  712 IF (ISYS.NE.0) GO TO 713
      WRITE (NTOUT,2080)
 2080 FORMAT (//'***', 35H INPUT ERROR: NO CELL CARD SUPPLIED)
      IERR=IERR+1
C
C     TEST INPUT ERROR COUNT
C
  713 IF (IERR.EQ.0) GO TO 715
      WRITE (NTOUT,2040) IERR
 2040 FORMAT (//' THERE ARE AT LEAST ',I4, 21H FATAL ERRORS IN THIS, 31H
     1 INPUT. EXECUTION IS TERMINATED)
      IERR=0
      IRNG=0
      GO TO 56
  715 IF (NRE.EQ.0) GO TO 720
      DO 324 L=1,NRE
      BEDING(L)=0.0
      KZ1=IZ1(L)
      DO 324 K=1,KZ1
      IF (K.EQ.LIND(L)) GO TO 324
      L1=MOD(KARI(L,K),3)
      IF (L1.EQ.0) L1=3
      L2=1+(KARI(L,K)-L1)/3
      KARI(L,K)=KVAR(L1,L2)
  324 CONTINUE
  716 IF (IRNG.EQ.0.AND.NC.EQ.0) GO TO 720
      CALL LINPAR
  720 IF (NAPID.EQ.0) GO TO 725
      IF (MCON.GT.0) GO TO 725
      WRITE (NTOUT,2065)
 2065 FORMAT (//'***NO CONNECTIVITY CARDS (E.G. TETCON CARDS)',
     1' SUPPLIED. NO APID RUN POSSIBLE')
      IZEIL=IZEIL+2
      NAPID=0
  725 RETURN
C***  FINISH THE PROGRAM: EXTRA <CR> REQUESTED FROM USER FOR PC'S  (RJA JAN02)
  611 WRITE (NTOUT,1021)
 1021 FORMAT (//'*** DLS 76 FINISHED ***'//)
	WRITE(6,1021)
	WRITE(6,'(''$ ***** HIT A <cr> TO CLOSE WINDOW'')')
	READ(5,'(1x)')
      REWIND NFILEA
      STOP
      END
      SUBROUTINE SYMOP
C     PURPOSE... DECODES THE SYMMETRY SPECIFICATIONS ON THE SYMEQ
C     OR ATOM CARDS
      DIMENSION IZEICH(18), IIP(6), ITEXT(2,2), XX(3)
      REAL*8 ID
      REAL IDEC,IZEICH
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK7/ IDEC(100),KA(2,150),LCELL(6)
      COMMON /BLOCK8/ SI(3,150),SIGN(3,3,150),SYMOPS(10,200),B(3,3,200),
     1BI(3,200),IND,MTR(200)
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
C BK     1IDIMM
      DATA ITEXT/'ATOM','SYME',' ','Q'/
      DATA IZEICH/'1','2','3','4','5','6','7','8','9','0','/','+','-',
     1            ',','X','Y','Z',' '/
 1000 FORMAT (//'***FORMAT ERROR IN SYMMETRY SPECIFICATIONS ON ',A4,
     1A2,' CARD NO ',I3,  1H:/1X,15A4/)
      KERR=IERR
      NNN=N
      IF (IND.EQ.2) NNN=NEQU
      DO 100 K=1,40
      I=42-K
      IF (IDEC(I).EQ.IZEICH(18)) GO TO 100
      IIP(6)=I
      GO TO 1
  100 CONTINUE
      IERR=IERR+1
    1 LL=1
      NH=IIP(6)
C     REMOVE BLANKS BETWEEN SYMBOLS
      DO 101 I=2,NH
      IF (IDEC(I).EQ.IZEICH(18)) GO TO 101
      LL=LL+1
      IDEC(LL)=IDEC(I)
  101 CONTINUE
      IDEC(LL+1)=IZEICH(18)
      IIP(6)=LL
      KK=LL+2
      DO 102 K=2,LL
      I=KK-K
      IF (IDEC(I).NE.IZEICH(14)) GO TO 102
      IIP(5)=I+1
      IIP(4)=IIP(5)-2
      GO TO 2
  102 CONTINUE
      IERR=IERR+1
    2 NH=IIP(4)
      KK=NH+2
      DO 103 K=2,NH
      I=KK-K
      IF (IDEC(I).NE.IZEICH(14)) GO TO 103
      IIP(3)=I+1
      IIP(2)=IIP(3)-2
      GO TO 3
  103 CONTINUE
      IERR=IERR+1
    3 IIP(1)=2
      LP1=-1
C     DECODE SYMBOLS
      DO 110 LP=1,3
      LP1=LP1+2
      LL1=IIP(LP1)
      LL2=IIP(LP1+1)
C     HOMOGENOUS PART OF SYMMETRY OPERATION
C     X, Y OR Z
      DO 106 I=1,3
      L1=I+14
      DO 105 L=LL1,LL2
      IF (IDEC(L).NE.IZEICH(L1)) GO TO 105
      K1=L-1
      K2=1
C     COEFFICIENTS OF X, Y OR Z
      DO 200 J=2,9
      IF (IDEC(K1).NE.IZEICH(J)) GO TO 200
      K2=J
      K1=L-2
      GO TO 5
  200 CONTINUE
      GO TO 5
  105 CONTINUE
      SIGN(LP,I,NNN)=0.0
      GO TO 106
    5 IF (IDEC(K1).EQ.IZEICH(12).OR.IDEC(K1).EQ.IZEICH(14).OR.IDEC(K1).E
     1Q.IZEICH(18)) GO TO 6
      IF (IDEC(K1).NE.IZEICH(13)) GO TO 7
      SIGN(LP,I,NNN)=-1.0*FLOAT(K2)
      GO TO 106
    6 SIGN(LP,I,NNN)=FLOAT(K2)
      GO TO 106
    7 IERR=IERR+1
  106 CONTINUE
C     TRANSLATIONS OR FIXED COORDINATES
      DO 107 I=LL1,LL2
      IF (IDEC(I).NE.IZEICH(11)) GO TO 107
      I1=I
      GO TO 8
  107 CONTINUE
C     AS AN INTEGER L2 OR -L2
      DO 126 I=LL1,LL2
      DO 120 IN=1,10
      IF (IDEC(I).NE.IZEICH(IN)) GO TO 120
      L2=IN
      IF (IN.EQ.10) L2=0
      K1=I-1
      K2=I+1
      DO 121 IX=15,17
      IF (IDEC(K2).EQ.IZEICH(IX)) GO TO 126
  121 CONTINUE
      GO TO 122
  120 CONTINUE
  126 CONTINUE
      GO TO 125
  122 IF (IDEC(K2).NE.IZEICH(12).AND.IDEC(K2).NE.IZEICH(13).AND.IDEC(K2)
     1.NE.IZEICH(14).AND.IDEC(K2).NE.IZEICH(18)) GO TO 124
      IF (IDEC(K1).EQ.IZEICH(12).OR.IDEC(K1).EQ.IZEICH(14).OR.IDEC(K1).E
     1Q.IZEICH(18)) GO TO 123
      IF (IDEC(K1).NE.IZEICH(13)) GO TO 124
      SI(LP,NNN)=-1.0*FLOAT(L2)
      GO TO 110
  123 SI(LP,NNN)=FLOAT(L2)
      GO TO 110
  125 SI(LP,NNN)=0.0
      GO TO 110
  124 IERR=IERR+1
      GO TO 110
C     AS A FRACTION L2/L1 OR -L2/L1
    8 DO 108 L=1,9
      IF (IDEC(I1+1).NE.IZEICH(L)) GO TO 108
      L1=L
      GO TO 9
  108 CONTINUE
      IERR=IERR+1
    9 DO 109 L=1,9
      IF (IDEC(I1-1).NE.IZEICH(L)) GO TO 109
      L2=L
      GO TO 10
  109 CONTINUE
      IERR=IERR+1
   10 NH=I1-2
      IF (IDEC(NH).EQ.IZEICH(12).OR.IDEC(NH).EQ.IZEICH(14).OR.IDEC(NH).E
     1Q.IZEICH(18)) GO TO 11
      IF (IDEC(NH).NE.IZEICH(13)) GO TO 12
      SI(LP,NNN)=-FLOAT(L2)/FLOAT(L1)
      GO TO 110
   11 SI(LP,NNN)=FLOAT(L2)/FLOAT(L1)
      GO TO 110
   12 IERR=IERR+1
  110 CONTINUE
      IF (IND.EQ.2) GO TO 201
C     ATOMS ON ATOM CARDS: MATRICES B AND VECTORS BI USED IN
C     SUBROUTINES DISDER AND NEWPAR
      DO 202 I=1,3
      BI(I,N)=SI(I,N)
      DO 202 K=1,3
  202 B(I,K,N)=SIGN(I,K,N)
C     CALCULATE REMAINING COORDINATES OF ATOMS IN SPECIAL POSITION
      DO 203 I=1,3
  203 XX(I)=X(I,N)
      DO 204 I=1,3
      XI=0.0
      DO 205 K=1,3
  205 XI=XI+B(I,K,N)*XX(K)
  204 X(I,N)=XI+BI(I,N)
      IF (KERR.NE.IERR) GO TO 208
      RETURN
  201 DETSIG=SIGN(1,1,NEQU)*SIGN(2,2,NEQU)*SIGN(3,3,NEQU)+SIGN(1,2,NEQU)
     1*SIGN(2,3,NEQU)*SIGN(3,1,NEQU)+SIGN(1,3,NEQU)*SIGN(2,1,NEQU)*SIGN(
     23,2,NEQU)-SIGN(1,2,NEQU)*SIGN(2,1,NEQU)*SIGN(3,3,NEQU)-SIGN(1,1,NE
     3QU)*SIGN(2,3,NEQU)*SIGN(3,2,NEQU)-SIGN(1,3,NEQU)*SIGN(2,2,NEQU)*SI
     4GN(3,1,NEQU)
      HN=DETSIG
      HN=ABS(HN)
      HN=ABS(HN-1.0)
      DD=0.1
      IF (HN.LT.DD) GO TO 13
      WRITE (NTOUT,1001) NEQU,DETSIG
 1001 FORMAT (//36H*** ERROR: DETERMINANT OF SYMEQ CARD,I3,3H IS,F8.3)
      IERR=IERR+1
      KERR=KERR+1
      GO TO 210
   13 K1=KA(1,NEQU)
C     COORDINATES OF SYMMETRICALLY DEPENDENT ATOMS
      DO 111 I=1,3
      X(I,N)=SIGN(I,1,NEQU)*X(1,K1)+SIGN(I,2,NEQU)*X(2,K1)+SIGN(I,3,NEQU
     1)*X(3,K1)+SI(I,NEQU)
  111 CONTINUE
C     ATOMS ON SYMEQ CARDS: MATRICES B
      DO 206 I=1,3
      DO 206 K=1,3
      BIK=0.0
      DO 207 L=1,3
  207 BIK=BIK+SIGN(I,L,NEQU)*B(L,K,K1)
  206 B(I,K,N)=BIK
  210 IF (KERR.NE.IERR) GO TO 208
      RETURN
  208 WRITE (NTOUT,1000) ITEXT(IND,1),ITEXT(IND,2),NNN,(SYMOPS(J,N),J=1,
     110)
      RETURN
      END
      SUBROUTINE DEFVAR
C      PURPOSE... COUNTS THE DIFFERENT KINDS OF VARIABLES
C              (COORDINATES,CELL,DISTANCES)
C         CALCULATES DIFFERENT ARRAYS TO IDENTIFY THE VARIABLES
      REAL*8 ID
      REAL IHR,IFAK,IGIN
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK9/ DELV(300),IFELD(6),IV(300),IV1(300),KVAR(3,200)
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
C BK     1IDIMM
      K1=0
C***   COORDINATES
      IF (N1.GT.0) GO TO 1
      WRITE (NTOUT,1000)
      IERR=IERR+1
      RETURN
    1 DO 100 I=1,N1
      DO 100 K=1,3
      IF (IFAK(K,I).NE.IHR) GO TO 2
      K1=K1+1
      KVAR(K,I)=K1
      IV(K1)=K
      IV1(K1)=I
      GO TO 100
    2 KVAR(K,I)=0
  100 CONTINUE
      NVO=K1
C***   CELL
      IF (IGIT.EQ.0) GO TO 3
      DO 101 I=1,6
      IF (IGIN(I).NE.IHR) GO TO 101
      K1=K1+1
      IV(K1)=I+3
      L1=K1-NVO
      IFELD(L1)=I
  101 CONTINUE
    3 NVG=K1-NVO
      NV1=NVG+NVO
      NV=NV1+NVD
      RETURN
 1000 FORMAT (//'***', 31H ERROR: THERE ARE NO ATOM CARDS)
      END
      SUBROUTINE DATOUT
C      PURPOSE... PRINTS THE CHECKED INPUT DATA
      INTEGER BT
	CHARACTER*90 STRING
      DIMENSION KXK(20), XKX(20), A2(6), LABEL(27), ISG(20)
      DIMENSION TXT(13), DTXT1(13), DTXT2(7), DTXT3(12), DTXT4(3), MX(3,
     12),DTXT5(28),DTXT6(4)
      DIMENSION IDXX(3), KXK1(20), IC(20)
      REAL*8 ID,MX,KXK,LABEL
      REAL IC,IBL,ICX,IHI,IHR,ISG,IDXX,IFAK,IGIN,ISTA,ISTB,KXK1,ISYMB
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK4/ RN,RSIGMA,RWERT,ICVTST
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK7/ IDEC(100),KA(2,150),LCELL(6)
      COMMON /BLOCK8/ SI(3,150),SIGN(3,3,150),SYMOPS(10,200),B(3,3,200),
     1BI(3,200),IND,MTR(200)
      COMMON /BLOCK9/ DELV(300),IFELD(6),IV(300),IV1(300),KVAR(3,200)
      COMMON S(50000)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
      COMMON /BLOC13/ ISYS,CVGTST,LSYS
      COMMON /BLOC14/ IZ1(35),WR(35),BEDING(35),RESKO(35,20),KARI(35,20)
     1,LINH(35),LIND(35)
      COMMON /BLOC15/ NAPID,IRNG,KRNG
      COMMON /BLOC16/ ITYPE(200),BT(15),IA,NZA,ICON(40,11),MD(40,18),NCA
     1PID,CAPID(15,4),IAT,IAB,IBT(40,4),DOT,WEIGHT(15,3),DELTAX(3),DELTA
     20(3),DA,DB,ANG,TOT(40,4),OTO(40,6),PI180,IMOV(3),MSYMAB(200),DOBIN
     3(400),KTYPE(20),IDK
C BK                                  ?????????
C BK      COMMON /BLOC17/ JS,MS,ITEST,IEQD(400),MEFF,WRF,IERR,MEQD,MCON,WFAC
C BK     1     ,IDIMM                 ?????????
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
      COMMON /BLOC18/ ROV,NC,NA,NS,IVF,IRNG1
      DATA IDXX,ISTA,ISTB,IBL/'DX (','DY (','DZ (','*','0 =',' '/
      DATA DTXT1/'NEW ','GAUS','S   ','  FU','LL  ','NO  ','  NO','    '
     1           ,' NO ','    ',' NO ','    ',' NO '/
      DATA DTXT2/'    ',' CAR','DS  ','YES ',' NO ','    ','NO  '/
      DATA DTXT3/'FIRS','T/LA','ST  ','FIRS','T/LA', 'ST  ','  NO',
     1     '    ','  NO','    ','  NO','    '/
      DATA MX/'  ','  ','  ',' X',' Y',' Z'/
      DATA DTXT4/'COEF','FICI','ENTS'/
      DATA DTXT5/'RAPH',' NG ','- NR','    ','DIAG','ONAL','SETS',' CYC'
     1          ,'LES ','NG /','  NR','CYCL','ES  ','YES ',' /  ','  / '
     2          ,' NO ','ALON','G   ','   M','ATRI','X   ','  EV','ERY '
     3          ,'FINA','L   ','EVER','Y   '/
      DATA DTXT6/'-','+','D','C'/
      IBED=0
      IF (IZEIL.EQ.0) GO TO 80
      CALL KOPF
   80 WRITE (NTOUT,999)
  999 FORMAT ( 23H PROGRAM SPECIFICATIONS/ 23H ----------------------/)
      WRITE (NTOUT,1000)
 1000 FORMAT (4X, 10HREFINEMENT,8X,  6HMATRIX,10X,  6HRANDOM,8X, 10HCONS
     1T DIST,5X, 28HCONST  RATIO   PRESCRIBED DO,8X,  4HCELL/7X,  4HTYPE
     2,9X,  9HINVERSION,7X, 11HCOORDINATES,5X,2( 10HREFINEMENT,6X), 10HA
     3DJUSTMENT,6X, 10HREFINEMENT//)
      DO 250 I=1,13
  250 TXT(I)=DTXT1(I)
      IF (INEW.EQ.0) GO TO 310
      IF (NRAGA.NE.0) GO TO 305
      TXT(2)=DTXT5(1)
      GO TO 310
  305 TXT(1)=DTXT5(2)
      TXT(2)=DTXT5(3)
      TXT(3)=DTXT5(4)
  310 IF (NDIAG.EQ.0) GO TO 315
      TXT(4)=DTXT5(5)
      TXT(5)=DTXT5(6)
  315 IF (IRNG.EQ.0) GO TO 320
      WRITE (NTOUT,1315) IRNG
 1315 FORMAT (36X,I2)
      TXT(6)=DTXT5(7)
  320 IF (NCYCLE.EQ.0) GO TO 330
      IF (NRAGA.NE.0) GO TO 325
      WRITE (NTOUT,1320) NCYCLE
 1320 FORMAT (51X,I2)
      TXT(7)=DTXT5(8)
      TXT(8)=DTXT5(9)
      GO TO 330
  325 NR=NCYCLE-NRAGA
      WRITE (NTOUT,1325) NRAGA,NR
 1325 FORMAT (51X,I2,4X,I2)
      TXT(7)=DTXT5(10)
      TXT(8)=DTXT5(11)
  330 IF (IDOB.EQ.0) GO TO 335
      WRITE (NTOUT,1330) IDOB
 1330 FORMAT (67X,I2)
      TXT(9)=DTXT5(12)
      TXT(10)=DTXT5(13)
  335 IF (NAPID.EQ.0) GO TO 340
      WRITE (NTOUT,1335) NAPID
 1335 FORMAT (83X,I2)
      TXT(11)=DTXT5(12)
      TXT(12)=DTXT5(13)
  340 IF (IGIT.EQ.0) GO TO 345
      TXT(13)=DTXT5(14)
  345 WRITE (NTOUT,1345) (TXT(I),I=1,13)
 1345 FORMAT (3X,3A4,5X,2A4,11X,A4,10X,2A4,9X,2A4,8X,2A4,8X,A4)
      WRITE (NTOUT,1400)
 1400 FORMAT (///' PARAMETERS FROM   CONVERGENCE',5X, 12HUSE VARIABLE,
     14X, 12HTRANSL FINAL/ 16H (CYC./APID/SET),6X,  4HTEST,9X, 12HDAMP
     2FACTOR,4X, 11HCOORDINATES//)
      DO 400 I=1,7
  400 TXT(I)=DTXT2(I)
      IF (NC+NA+NS.EQ.0) GO TO 405
      WRITE (NTOUT,1405) NC,NA,NS
 1405 FORMAT (1X,I4,2I5)
      TXT(2)=DTXT5(15)
      TXT(3)=DTXT5(16)
  405 IF (ICAL.EQ.0) GO TO 410
      TXT(4)=DTXT5(17)
  410 IF (IVF.EQ.0) GO TO 415
      TXT(5)=DTXT5(14)
  415 IF (IMOV(1)+IMOV(2)+IMOV(3).EQ.0) GO TO 420
      WRITE (NTOUT,1410) (MX(I,IMOV(I)+1),I=1,3)
 1410 FORMAT (56X,3A2)
      TXT(6)=DTXT5(18)
      TXT(7)=DTXT5(19)
  420 WRITE (NTOUT,1450) (TXT(I),I=1,7)
 1450 FORMAT (1X,3A4,9X,A4,13X,A4,9X,2A4)
      WRITE (NTOUT,1499)
 1499 FORMAT (/// ' OUTPUT CONTROL'/ 15H --------------)
      DO 425 I=1,3
  425 TXT(I)=DTXT4(I)
      IF (ICOR.GE.0) GO TO 430
      TXT(1)=DTXT5(20)
      TXT(2)=DTXT5(21)
      TXT(3)=DTXT5(22)
  430 WRITE (NTOUT,1500) (TXT(I),I=1,3),NFILEA
 1500 FORMAT (  //6X,  4HLIST,12X,  4HLIST,8X, 11HLIST MATRIX,5X, 11HL
     1IST CORREL,4X, 16HWRITE PARAMETERS/4X, 10HPARAMETERS,7X,  9HDISTAN
     2CES,5X, 10HAND VECTOR,6X,3A4,3X, 11HON NFILEA (,I2,  1H)/)
      DO 500 I=1,12
  500 TXT(I)=DTXT3(I)
      IF (IORT.EQ.0) GO TO 505
      TXT(1)=DTXT5(23)
      TXT(2)=DTXT5(24)
      TXT(3)=DTXT5(4)
  505 IF (IDIS.EQ.0) GO TO 510
      TXT(4)=DTXT5(23)
      TXT(5)=DTXT5(24)
      TXT(6)=DTXT5(4)
  510 IF (IMAT-1) 525,515,520
  515 TXT(7)=DTXT5(25)
      TXT(8)=DTXT5(26)
      GO TO 525
  520 TXT(7)=DTXT5(27)
      TXT(8)=DTXT5(28)
  525 IF (IABS(ICOR)-1) 540,530,535
  530 TXT(9)=DTXT5(25)
      TXT(10)=DTXT5(26)
      GO TO 540
  535 TXT(9)=DTXT5(27)
      TXT(10)=DTXT5(28)
  540 IF (IPCH-1) 555,545,550
  545 TXT(11)=DTXT5(25)
      TXT(12)=DTXT5(26)
      GO TO 555
  550 TXT(11)=DTXT5(27)
      TXT(12)=DTXT5(28)
  555 WRITE (NTOUT,1580) (TXT(I),I=1,12)
 1580 FORMAT (  1H ,3X,3A4,4X,3A4,6X,2A4,8X,2A4,8X,2A4)
      WRITE (NTOUT,1600)
 1600 FORMAT (/////' CONTROL FACTORS'/ 16H ---------------//9X, 15HDAM
     1PING FACTORS,11X, 12HSTARTING  NR,4X, 59HCONVERGENCE    LINRES WEI
     2GHTS  WTS FACTOR FOR   CORR MATRIX/ 62H COORDINATES   CELL    DIST
     3ANCES   RANDOM GENER    TEST FACTOR,4X, 44HMULTIPLIED BY  OUT OF R
     4ANGE D'S  TEST FACTOR/)
      SQW=SQRT(WFAC)
      CORRL=CORR
      IF (ICOR.LT.0) CORRL=0.0
      WRITE (NTOUT,1605) (FUDGE(I),I=1,3),IRNG1,CVGTST,WRF,SQW,CORRL
 1605 FORMAT (1X,F8.3,F10.3,F11.3,I13,F18.7,F18.5,F13.2,F14.2////)
      NRS=NRE-NRH
      WRITE (NTOUT,1017) N1,N2,M,NV,NRH,NRS
 1017 FORMAT ( 23H PRELIMINARY STATISTICS/ 23H ----------------------//2
     1(3X,  9HNUMBER OF,5X),3(1X,  9HNUMBER OF,6X)/ 15H INDEPEND.ATOMS,4
     2X, 12HDEPEND.ATOMS,4X,  9HDISTANCES,7X,  9HVARIABLES,6X, 12HRESTRI
     3CTIONS//2(5X,I4,8X),3X,I4,11X,I4,8X,I2,  6H HARD/,I2,  5H SOFT)
      CALL KOPF
      GO TO 90
      ENTRY PAREX
      IBED=1
      IF (IZEIL.NE.0) CALL KOPF
      WRITE (NTOUT,1001)
 1001 FORMAT ( 17H FINAL PARAMETERS/ 17H ----------------/)
   90 DO 40 K=1,3
   40 A2(K)=A(K)
      DO 45 K=4,6
   45 A2(K)=A(K)/PI180
      IF (NVG.EQ.0) GO TO 50
      K1=NVO+1
      K2=NV1
      GO TO 51
   50 K1=0
      K2=0
   51 IF (IBED.EQ.1) GO TO 52
      WRITE (NTOUT,1002)
   52 WRITE (NTOUT,1006) (A2(I),IGIN(I),I=1,6),K1,K2
 1002 FORMAT ( 19H INITIAL PARAMETERS/ 19H ------------------//20X,  8HR
     1=REFINE,6X, 33HD=DEPENDENT IN A SPECIAL POSITION//20X, 11HI=INVARI
     2ANT,3X, 39HL=DEPENDENT DUE TO A LINEAR RESTRICTION/)
 1006 FORMAT (1X,15HCELL PARAMETERS///6X,1HA,11X,1HB,11X,1HC,10X,5HALPHA
     1,8X,4HBETA,7X,5HGAMMA,7X,9HVARIABLES//3(F10.4,1X,A1),3(F10.2,1X,A1
     2),5X,I3,2H -,I3///,1X,15HATOM PARAMETERS//1X,4HATOM,6X,1HX,9X,1H
     3Y,9X,1HZ,9X,4HTYPE,6X,1HN,6X,17HSPECIAL POSITIONS/)
      IZEIL=23
      IF (IBED.EQ.1) IZEIL=19
      DO 102 I=1,N1
      CALL KOPF1
      DO 101 K=1,3
      IF (KVAR(K,I).EQ.0) GO TO 101
      K1=KVAR(K,I)
      GO TO 105
  101 CONTINUE
      K1=0
  105 STRING=' '		!CLEAR IT
	WRITE (STRING,1007) ID(I),(X(K,I),IFAK(K,I),K=1,3),KTYPE(ITYPE(I)),
     1K1
 1007 FORMAT (1X,A6,3(F8.5,1X,A1),5X,A3,4X,I3)

C  105 WRITE (NTOUT,1007) ID(I),(X(K,I),IFAK(K,I),K=1,3),KTYPE(ITYPE(I)),
C     1K1
C 1007 FORMAT (1X,A6,3(F8.5,1X,A1),5X,A3,4X,I3)
      IF (MTR(I).NE.0)THEN
		WRITE (STRING(59:89),1050) (SYMOPS(K,I),K=1,8)
1050		FORMAT(7A4,A2)
	ENDIF
	WRITE(NTOUT,'(A90)')STRING
C      WRITE (NTOUT,1050) (SYMOPS(K,I),K=1,8)
C 1050 FORMAT (1H+,58X,7A4,A2)  !RJA 

  102 CONTINUE
      IF (IBED.EQ.1) GO TO 99
      WRITE (NTOUT,1005)
      IZEIL=IZEIL+2
 1005 FORMAT ( ' (N=NO OF THE FIRST VAR. IN THIS LINE)')
   99 WRITE (NTOUT,1004)
      IZEIL=IZEIL+2
      CALL KOPF1
 1004 FORMAT (//)
      IF (N.EQ.N1) GO TO 34
 1008 FORMAT (1H )
      CALL KOPF1
      WRITE (NTOUT,1009)
 1009 FORMAT (1X,29HSYMMETRICALLY DEPENDENT ATOMS,18X,19HSYMMETRY OPERAT
     1IONS)
      CALL KOPF1
      WRITE (NTOUT,1008)
      NH=N1+1
      DO 104 I=NH,N
      CALL KOPF1
  104 WRITE (NTOUT,1010) ID(I),(X(K,I),K=1,3),(SYMOPS(K,I),K=1,10)
 1010 FORMAT (1X,A6,3(F8.5,2X),11X,10A4)
   34 IF (IBED.EQ.1) RETURN
      WRITE (NTOUT,1008)
      CALL KOPF1
      IF (NRE.EQ.0) GO TO 200
      IF (IZEIL+7.GT.LETZTE) CALL KOPF
      WRITE (NTOUT,1014)
 1014 FORMAT (' LINEAR RESTRICTIONS'/  ' NR',5X,  6HWEIGHT/)
      IZEIL=IZEIL+4
      CALL KOPF1
      DO 108 L=1,NRE
      L1=IZ1(L)
      DO 107 I=1,L1
      IF (LINH(L).EQ.0.OR.I.NE.LIND(L)) GO TO 300
      LAT2=MOD(KARI(L,I),3)
      IF (LAT2.EQ.0) LAT2=3
      LAT1=1+(KARI(L,I)-LAT2)/3
      GO TO 301
  300 K1=KARI(L,I)
      LAT1=IV1(K1)
      LAT2=IV(K1)
  301 KXK(I)=ID(LAT1)
      KXK1(I)=IDXX(LAT2)
      ISG(I)=IBL
      IF (RESKO(L,I).LT.-1.E-8) ISG(I)=DTXT6(1)
      IF (RESKO(L,I).GT.1.E-8.AND.I.NE.1) ISG(I)=DTXT6(2)
  107 XKX(I)=ABS(RESKO(L,I))
      IF (LINH(L).EQ.0) GO TO 302
      WRITE (NTOUT,1085) L,ISTB,(IBL,ISG(K),XKX(K),ISTA,KXK1(K),KXK(K),K
     1=1,L1)
 1085 FORMAT (1X,I2,6X,  4HHARD,6X,A3,A1,4(A1,F4.1,A1,A4,A6,  1H),A1),3(
     1/24X,4(A1,F4.1,A1,A4,A6,  1H),A1)))
      GO TO 303
  302 WRQ=SQRT(WR(L))/WRF
      WRITE (NTOUT,1015) L,WRQ,ISTB,(IBL,ISG(K),XKX(K),ISTA,KXK1(K),KXK(
     1K),K=1,L1)
  303 IF (L1.GT.4) IZEIL=IZEIL+1
      CALL KOPF1
 1015 FORMAT (1X,I2,E12.3,4X,A3,A1,4(A1,F4.1,A1,A4,A6,  1H),A1),3(/24X,4
     1(A1,F4.1,A1,A4,A6,  1H),A1)))
  108 CONTINUE
      WRITE (NTOUT,1008)
      CALL KOPF1
  200 IF (NFUNC.EQ.0) GO TO 115
C     PRINT OUT BONDIS FUNCTIONS
      IF (IZEIL+6+2*IA.GT.LETZTE) CALL KOPF
      WRITE (NTOUT,1021)
 1021 FORMAT (' BOND DISTANCE FUNCTIONS'//'  BOND - TYPE',37X,  8HFU
     1NCTION,39X,  6HW(T-O),4X,  6HW(O-O),4X,  6HW(T-T)/)
      IZEIL=IZEIL+5
      CALL KOPF1
      DO 110 I=1,IA
      IB3=BT(I)/(IDK*IDK)
      IB1=BT(I)-IB3*IDK*IDK
      IB2=IB1/IDK
      IB1=IB1-IB2*IDK
      IB1=KTYPE(IB1)
      IB2=KTYPE(IB2-2)
      IB3=KTYPE(IB3-2)
      WRITE (NTOUT,1023) IB1,IB2,IB3,CAPID(I,1),CAPID(I,2),CAPID(I,4),CA
     1PID(I,3),CAPID(I,4),(WEIGHT(I,L),L=1,3)
 1023 FORMAT (1X,A4,  2H- ,A3,  4H... ,A4,  5H DO =,F8.5,  2H +,E12.4,9H
     1 * (TOT -,F6.1,  3H) +,E12.4,9H * (TOT -,F6.1,4H)**2,3F10.4)
  110 CALL KOPF1
      WRITE (NTOUT,1008)
      CALL KOPF1
      DO 210 I=1,IA
      IB3=BT(I)/(IDK*IDK)
      IB1=BT(I)-IB3*IDK*IDK
      IB2=IB1/IDK
      IB1=IB1-IB2*IDK
      IB1=KTYPE(IB1)
      IB3=KTYPE(IB3-2)
      WRITE (NTOUT,1040) IB1,IB3,CAPID(I,4)
 1040 FORMAT (  5H DO( ,A3,  2H- ,A3, 33H) IS CALCULATED USING AN ANGLE
     1OF,F6.1,  8H DEGREES)
  210 CALL KOPF1
      WRITE (NTOUT,1008)
      CALL KOPF1
  115 IF (IDOB.EQ.0) GO TO 120
      IF (IZEIL+6+NVD.GT.LETZTE) CALL KOPF
      WRITE (NTOUT,1011)
 1011 FORMAT (/' PRESCRIBED DISTANCE VARIABLES'/  //3X,  2HDO,8X,  1
     1HN,5X, 28H(D=CONSTANT RATIO TO VAR. N)/)
      IZEIL=IZEIL+6
      DO 100 K=1,NDO
      IF (LDR(K).EQ.0) GO TO 100
      IC(LDR(K))=IHR
  100 CONTINUE
      DO 106 K=1,NDO
      CALL KOPF1
      IF (LDR(K).EQ.0) GO TO 33
      K1=NV1+LDR(K)
      ICX=IC(LDR(K))
      IC(LDR(K))=DTXT6(3)
      GO TO 106
   33 K1=0
      ICX=IHI
  106 WRITE (NTOUT,1013) DOB(K),ICX,K1
 1013 FORMAT (F7.4,1X,A1,I6)
  120 IF (MCON.EQ.0) GO TO 160
      WRITE (NTOUT,1008)
      IZEIL=IZEIL+1
      IF (IZEIL+20.GT.LETZTE) CALL KOPF
      IZEIL=IZEIL+2
      WRITE (NTOUT,1060)
 1060 FORMAT ( 35H CONNECTIVITY OF TETRAHEADRAL ATOMS/)
      DO 150 NZ1=1,NZA,3
      NZ2=NZ1+2
      IF (NZ2.GT.NZA) NZ2=NZA
      IF (IZEIL+18.GT.LETZTE) CALL KOPF
      L1=-3
      DO 145 I=1,9
      L1=L1+3
      L=L1
      DO 145 NZ=NZ1,NZ2
      L=L+1
      ICO=ICON(NZ,I)
      IF (ICO.EQ.0) GO TO 140
      IF (ICO.LT.0) ICO=-ICO
      LABEL(L)=ID(ICO)
      GO TO 145
  140 LABEL(L)=MX(1,1)
  145 CONTINUE
      IT=NZ2-NZ1+1
      WRITE (NTOUT,1065) (LABEL(18+I),I=1,IT)
 1065 FORMAT (  //16X,A6,37X,A6,37X,A6)
      WRITE (NTOUT,1065) (IHI,I=1,IT)
      WRITE (NTOUT,1065) (LABEL(6+I),I=1,IT)
      WRITE (NTOUT,1065) (IHI,I=1,IT)
      WRITE (NTOUT,1075) (LABEL(15+I),LABEL(3+I),LABEL(I),LABEL(9+I),LAB
     1EL(21+I),I=1,IT)
 1075 FORMAT (  //4(A6,  2H- ),A6,5X,4(A6,  2H- ),A6,5X,4(A6,  2H- ),A
     16)
      WRITE (NTOUT,1065) (IHI,I=1,IT)
      WRITE (NTOUT,1065) (LABEL(12+I),I=1,IT)
      WRITE (NTOUT,1065) (IHI,I=1,IT)
      WRITE (NTOUT,1065) (LABEL(24+I),I=1,IT)
      WRITE (NTOUT,1008)
      IZEIL=IZEIL+17
      CALL KOPF1
  150 CONTINUE
  160 IF (MEQD.EQ.0) GO TO 165
      WRITE (NTOUT,1030)
 1030 FORMAT (//' THE FOLLOWING DISTANCES, SPECIFIED ON A CONNECTIVITY'
     1, 44H CARD HAVE BEEN REDEFINED WITH A DISTAN CARD/ '  DISTANCE',
     211X, 26HNEW DO AND WEIGHT(SQUARED)/)
      IZEIL=IZEIL+3
      DO 203 J=1,M
      IF (ISYMB(J).NE.DTXT6(4)) GO TO 203
      WRITE (NTOUT,1031) ID(NN(J,1)),ID(NN(J,2)),DOB(J),W(J),ISYMB(J)
      CALL KOPF1
 1031 FORMAT (1X,A6,  2H- ,A6,2F15.4,1X,A1)
  203 CONTINUE
      IF (NAPID.EQ.0) GO TO 165
      WRITE (NTOUT,1035)
 1035 FORMAT (//' DISTANCES MARKED WITH C WILL NOT BE CHANGED', 17H IN A
     1N APID CYCLE//)
      IZEIL=IZEIL+3
      CALL KOPF1
  165 IF (NV.LE.M) GO TO 170
      WRITE (NTOUT,1080)
 1080 FORMAT (//'***', 47H ERROR: THERE ARE MORE VARIABLES THAN EQUATIONS
     1/  //11X, 35HEXECUTION OF THIS JOB IS TERMINATED)
      IERR=IERR+1
  170 RETURN
      END
      SUBROUTINE SPEICH
C      PURPOSE... STORES SINUS AND COSINUS OF ALPHA,BETA,GAMMA
C            CALCULATES THE METRIC TENSOR AND ITS FIRST AND SECOND
C            DERIVATIVES.
      REAL*8 ID
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON S(50000)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
C BK      COMMON /BLOC13/ ISYS,CVGTST
      COMMON /BLOC13/ ISYS,CVGTST,LSYS
      DO 100 K=1,3
      DO 100 L=1,6
      GABL1(1,K,L)=0.0
      GABL1(2,K,L)=0.0
      GABL1(3,K,L)=0.0
      DO 100 LS=L,6
      GABL2A(K,L,LS)=0.0
      GABL2B(K,L,LS)=0.0
      GABL2C(K,L,LS)=0.0
  100 CONTINUE
C***  STORE SIN AND COS
      SIN4=SIN(A(4))
      SIN5=SIN(A(5))
      SIN6=SIN(A(6))
      COS4=COS(A(4))
      COS5=COS(A(5))
      COS6=COS(A(6))
C***  CALCULATE METRIC TENSOR G
      G(1,1)=A(1)**2
      G(2,2)=A(2)**2
      G(3,3)=A(3)**2
      G(1,2)=A(1)*A(2)*COS6
      G(1,3)=A(1)*A(3)*COS5
      G(2,3)=A(2)*A(3)*COS4
      IF (IGIT.EQ.0) GO TO 10
C***  FIRST DERIVATIVE OF G
      GO TO (1,2,3,4), ISYS
C***   TRICLINIC,MONOCLINIC,ORTHORHOMBIC
    1 GABL1(1,1,1)=2.0*A(1)
      GABL1(1,2,1)=A(2)*COS6
      GABL1(1,3,1)=A(3)*COS5
      GABL1(1,2,2)=A(1)*COS6
      GABL1(2,2,2)=2.0*A(2)
      GABL1(2,3,2)=A(3)*COS4
      GABL1(1,3,3)=A(1)*COS5
      GABL1(2,3,3)=A(2)*COS4
      GABL1(3,3,3)=2.0*A(3)
      GABL1(2,3,4)=-A(2)*A(3)*SIN4
      GABL1(1,3,5)=-A(1)*A(3)*SIN5
      GABL1(1,2,6)=-A(1)*A(2)*SIN6
      GO TO 5
C***   TRIGONAL,HEXAGONAL,TETRAGONAL
    2 GABL1(1,1,1)=2.0*A(1)
      GABL1(2,2,1)=GABL1(1,1,1)
      GABL1(1,2,1)=GABL1(1,1,1)*COS6
      GABL1(3,3,3)=2.0*A(3)
      GO TO 5
C***   RHOMBOHEDRAL
    3 GABL1(1,1,1)=2.0*A(1)
      GABL1(2,2,1)=GABL1(1,1,1)
      GABL1(3,3,1)=GABL1(1,1,1)
      GABL1(1,2,1)=GABL1(1,1,1)*COS4
      GABL1(1,3,1)=GABL1(1,1,1)*COS4
      GABL1(2,3,1)=GABL1(1,1,1)*COS4
      GABL1(1,2,4)=-A(1)*A(1)*SIN4
      GABL1(1,3,4)=GABL1(1,2,4)
      GABL1(2,3,4)=GABL1(1,2,4)
      GO TO 5
C***   CUBIC
    4 GABL1(1,1,1)=2.0*A(1)
      GABL1(2,2,1)=GABL1(1,1,1)
      GABL1(3,3,1)=GABL1(1,1,1)
C***  SECOND DERIVATIVE OF G
    5 IF (INEW.EQ.0) GO TO 10
      GO TO (6,7,8,9), ISYS
C***   TRICLINIC,MONOCLINIC,ORTHORHOMBIC
    6 GABL2A(1,1,1)=2.0
      GABL2A(2,1,2)=COS6
      GABL2A(3,1,3)=COS5
      GABL2A(3,1,5)=-A(3)*SIN5
      GABL2A(2,1,6)=-A(2)*SIN6
      GABL2B(2,2,2)=2.0
      GABL2B(3,2,3)=COS4
      GABL2B(3,2,4)=-A(3)*SIN4
      GABL2A(2,2,6)=-A(1)*SIN6
      GABL2C(3,3,3)=2.0
      GABL2B(3,3,4)=-A(2)*SIN4
      GABL2A(3,3,5)=-A(1)*SIN5
      GABL2B(3,4,4)=-G(2,3)
      GABL2A(3,5,5)=-G(1,3)
      GABL2A(2,6,6)=-G(1,2)
      GO TO 10
C***   TRIGONAL,HEXAGONAL,TETRAGONAL
    7 GABL2A(1,1,1)=2.0
      GABL2B(2,1,1)=2.0
      GABL2A(2,1,1)=2.0*COS6
      GABL2C(3,3,3)=2.0
      GO TO 10
C***   RHOMBOHEDRAL
    8 GABL2A(1,1,1)=2.0
      GABL2B(2,1,1)=2.0
      GABL2C(3,1,1)=2.0
      GABL2A(2,1,1)=2.0*COS4
      GABL2A(3,1,1)=GABL2A(2,1,1)
      GABL2B(3,1,1)=GABL2A(2,1,1)
      GABL2A(2,4,4)=-A(1)*A(1)*COS4
      GABL2A(3,4,4)=GABL2A(2,4,4)
      GABL2B(3,4,4)=GABL2A(2,4,4)
      GABL2A(2,1,4)=-GABL1(1,1,1)*SIN4
      GABL2A(3,1,4)=GABL2A(2,1,4)
      GABL2B(3,1,4)=GABL2A(2,1,4)
      GO TO 10
C***   CUBIC
    9 GABL2A(1,1,1)=2.0
      GABL2B(2,1,1)=2.0
      GABL2C(3,1,1)=2.0
   10 DO 104 I=1,3
      DO 104 K=I,3
  104 G(K,I)=G(I,K)
      IF (IGIT.EQ.0) RETURN
      DO 101 I=1,3
      DO 101 K=I,3
      DO 101 L=1,6
  101 GABL1(K,I,L)=GABL1(I,K,L)
      IF (INEW.EQ.0) RETURN
      DO 110 L=1,6
      DO 110 LS=L,6
      GABL2A(1,LS,L)=GABL2A(1,L,LS)
      GABL2B(1,L,LS)=GABL2A(2,L,LS)
      GABL2A(2,LS,L)=GABL2A(2,L,LS)
  110 GABL2B(1,LS,L)=GABL2A(2,LS,L)
      DO 111 L=1,6
      DO 111 LS=L,6
      GABL2C(1,L,LS)=GABL2A(3,L,LS)
      GABL2A(3,LS,L)=GABL2A(3,L,LS)
      GABL2C(1,LS,L)=GABL2A(3,LS,L)
  111 GABL2B(2,LS,L)=GABL2B(2,L,LS)
      DO 113 L=1,6
      DO 113 LS=L,6
      GABL2C(2,L,LS)=GABL2B(3,L,LS)
      GABL2B(3,LS,L)=GABL2B(3,L,LS)
      GABL2C(2,LS,L)=GABL2B(3,LS,L)
  113 GABL2C(3,LS,L)=GABL2C(3,L,LS)
      RETURN
      END
      SUBROUTINE DISDER
C      PURPOSE... CALCULATES DISTANCES AND THEIR DERIVATIVES,
C      THE NORMAL MATRIX AND VECTOR
      DIMENSION PDELX(3,3), ABL1(300), DELTAX(3), IMG(96), ALIN(6,90)
      DIMENSION GMALDX(3), GIKDAL(3,6), LRE(6), A1(3)
      REAL*8 ID
      REAL IFAK
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK4/ RN,RSIGMA,RWERT,ICVTST
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK8/ SI(3,150),SIGN(3,3,150),SYMOPS(10,200),B(3,3,200),
     1BI(3,200),IND,MTR(200)
      COMMON /BLOCK9/ DELV(300),IFELD(6),IV(300),IV1(300),KVAR(3,200)
      COMMON S(50000)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
      COMMON /BLOC12/ SYMSIG(3,3,40),KATOR(400),NN11(400,2)
      COMMON /BLOC14/ IZ1(35),WR(35),BEDING(35),RESKO(35,20),KARI(35,20)
     1,LINH(35),LIND(35)
      COMMON /BLOC15/ NAPID,IRNG,KRNG
C BK      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
      DATA CHAR/'L'/
      INDEX(I,K)=NV*I+K-(I*(I-1))/2
      RO=0.0
      DD=1.0E-15
C***  CLEAR NORMAL MATRIX AND VECTOR
      DO 97 I=1,NV
   97 T(I)=0.0
      L=(NV*(NV+3))/2
      DO 96 I=1,L
   96 S(I)=0.0
      IF (NRAGA.EQ.0) GO TO 51
      IF (NCYC.LE.NRAGA) GO TO 50
      INEW=1
      GO TO 51
   50 INEW=0
   51 CONTINUE
      DO 200 J=1,M
      BIL=2.0*W(J)
      NIA=NN(J,1)
      NHA=NN(J,2)
      DO 99 K=1,3
   99 DELTAX(K)=X(K,NIA)-X(K,NHA)
      DO 100 I=1,3
      DO 100 K=I,3
      PDELX(I,K)=DELTAX(I)*DELTAX(K)
      PDELX(K,I)=PDELX(I,K)
  100 CONTINUE
      IF (IGIT.EQ.0) GO TO 16
      DO 104 L=1,6
      DO 104 I=1,3
      SUM=0.0
      DO 103 K=1,3
  103 SUM=SUM+GABL1(I,K,L)*DELTAX(K)
  104 GIKDAL(I,L)=SUM
C***  CALCULATE DISTANCES
   16 SUM=0.0
      DO 102 I=1,3
      SUMM=0.0
      DO 101 K=1,3
      SUMM=SUMM+G(I,K)*DELTAX(K)
  101 SUM=SUM+G(I,K)*PDELX(I,K)
  102 GMALDX(I)=SUMM
      IF (SUM.GE.DD) GO TO 1
      WRITE (NTOUT,1001) J,SUM
 1001 FORMAT (//17H***WARNING: D**2(,I4,2H)=,F10.4)
    1 D(J)=SQRT(SUM)
      JJ=LJ(J)
      DELTAD(J)=DOB(JJ)-D(J)
      RO=RO+W(J)*DELTAD(J)*DELTAD(J)
      KK=0
      IF (NVO.EQ.0) GO TO 2
      NI=NN11(J,1)
      NH=NN11(J,2)
C     LRE(K) = L OR LRE(K+3) = L: COORDINATE K OF ATOM NI OR NH
C     DEPENDENT BY HARD LINEAR RESTRICTION L
      DO 302 K=1,3
      LRE(K)=0
      LRE(K+3)=0
      DO 302 L=1,NRE
      IF (LINH(L).EQ.0) GO TO 302
      IF (IFAK(K,NI).NE.CHAR) GO TO 303
      IF (KARI(L,LIND(L)).EQ.3*(NI-1)+K) LRE(K)=L
  303 IF (IFAK(K,NH).NE.CHAR.OR.KATOR(J).GT.0) GO TO 302
      IF (KARI(L,LIND(L)).EQ.3*(NH-1)+K) LRE(K+3)=L
  302 CONTINUE
C     STORE NUMBERS OF INDEPENDENT VARIABLES IN IMG(KK)
C     FROM RESTRICTIONS
      DO 304 K=1,3
      IF (LRE(K).EQ.0) GO TO 305
      IZ=IZ1(LRE(K))-1
      DO 306 I=1,IZ
      IF (I.EQ.LIND(LRE(K))) GO TO 306
      KK=KK+1
      IMG(KK)=KARI(LRE(K),I)
  306 CONTINUE
  305 IF (LRE(K+3).EQ.0.OR.KATOR(J).GT.0) GO TO 307
      IZ=IZ1(LRE(K+3))-1
      DO 308 I=1,IZ
      IF (I.EQ.LIND(LRE(K+3))) GO TO 308
      KK=KK+1
      IMG(KK)=KARI(LRE(K+3),I)
  308 CONTINUE
C     FROM NORMALLY REFINED COORDINATES
  307 IF (KVAR(K,NI).EQ.0) GO TO 309
      KK=KK+1
      IMG(KK)=KVAR(K,NI)
  309 IF (KVAR(K,NH).EQ.0.OR.KATOR(J).GT.0) GO TO 304
      KK=KK+1
      IMG(KK)=KVAR(K,NH)
  304 CONTINUE
C     ELIMINATE NUMBERS IN IMG THAT APPEAR MORE THAN ONCE
      KKR=KK-1
      DO 310 K=2,KK
  313 DO 311 KL=2,K
      IF (IMG(KL-1).NE.IMG(K)) GO TO 311
      IF (K.LE.KKR) GO TO 341
      KKR=KKR-1
      GO TO 315
  341 DO 312 KH=K,KKR
  312 IMG(KH)=IMG(KH+1)
      KKR=KKR-1
      GO TO 313
  311 CONTINUE
      IF (K.EQ.KKR+1) GO TO 315
  310 CONTINUE
  315 KK=KKR+1
      IM=6
      IF (KATOR(J).GT.0) IM=3
C     CALCULATE MATRIX ALIN(I,K) CONTAINING COEFFICIENTS OF THE
C     VARIABLES IMG(K) (INTERMEDIATE FORM OF ALIN(I,K))
      DO 316 I=1,IM
      IF (LRE(I).EQ.0) GO TO 317
C     CASE OF LINEAR RESTRICTIONS
      DO 318 K=1,KK
      ALIN(I,K)=0.0
      IZ=IZ1(LRE(I))
      DO 319 IL=1,IZ
      IF (IL.EQ.LIND(LRE(I))) GO TO 319
      IF (KARI(LRE(I),IL).NE.IMG(K)) GO TO 319
      IRL=IL
      GO TO 320
  319 CONTINUE
      GO TO 318
  320 ALIN(I,K)=-RESKO(LRE(I),IRL)/RESKO(LRE(I),LIND(LRE(I)))
  318 CONTINUE
      GO TO 316
C     CASE OF NORMALLY REFINED VARIABLES
  317 NI1=NI
      KI=I
      IF (I.LE.3) GO TO 321
      NI1=NH
      KI=I-3
  321 DO 301 K=1,KK
      ALIN(I,K)=0.0
      IF (KVAR(KI,NI1).NE.IMG(K)) GO TO 301
      ALIN(I,K)=1.0
  301 CONTINUE
  316 CONTINUE
C     CALCULATE THE COEFFICIENT MATRIX USED IN THE DERIVATIVES
      IF (KATOR(J).NE.0) GO TO 322
C     ATOM 1 AND ATOM 2 INDEPENDENT AND IN GENERAL POSITION
      DO 323 I=1,3
      DO 323 K=1,KK
  323 ALIN(I,K)=ALIN(I,K)-ALIN(I+3,K)
      GO TO 329
  322 IF (KATOR(J).LT.0) GO TO 324
C     ATOM1 AND ATOM2 ARE SYMMETRICALLY EQUIVALENT
      DO 325 K=1,KK
      DO 326 I=1,3
      A1(I)=0.0
      DO 326 L=1,3
  326 A1(I)=A1(I)+SYMSIG(I,L,KATOR(J))*ALIN(L,K)
      DO 325 I=1,3
  325 ALIN(I,K)=A1(I)
      GO TO 329
C     GENERAL CASE
  324 DO 327 K=1,KK
      DO 328 I=1,3
      A1(I)=0.0
      DO 328 L=1,3
  328 A1(I)=A1(I)+B(I,L,NIA)*ALIN(L,K)-B(I,L,NHA)*ALIN(L+3,K)
      DO 327 I=1,3
  327 ALIN(I,K)=A1(I)
C***  FIRST DERIVATIVESW.R.T. COORDINATES
  329 DO 330 K=1,KK
      L1=IMG(K)
      SUM=0.0
      DO 331 I=1,3
  331 SUM=SUM+ALIN(I,K)*GMALDX(I)
  330 ABL1(L1)=SUM/D(J)
    2 KMG=KK
C***  FIRST DERIVATIVES W.R.T. CELL PAR.
      IF (IGIT.NE.1) GO TO 3
      DO 106 L1=1,NVG
      LH=IFELD(L1)
      L=L1+NVO
      SUM=0.0
      DO 105 I=1,3
      DO 105 K=1,3
  105 SUM=SUM+PDELX(I,K)*GABL1(I,K,LH)
      ABL1(L)=0.5*SUM/D(J)
      KMG=KMG+1
      IMG(KMG)=L
  106 CONTINUE
C***  CALCULATE VECTOR
    3 DO 107 KA=1,KMG
      I=IMG(KA)
      T(I)=T(I)+BIL*DELTAD(J)*ABL1(I)
  107 CONTINUE
      IF (INEW.EQ.1) GO TO 4
C***  CALCULATE NEWTON-GAUSS NORMAL MATRIX
      DO 110 KA=1,KMG
      L1=IMG(KA)
      DO 110 KB=KA,KMG
      L2=IMG(KB)
      IF (L1.GT.L2) GO TO 108
      KLM=INDEX(L1,L2)
      GO TO 109
  108 KLM=INDEX(L2,L1)
  109 PABL=ABL1(L1)*ABL1(L2)
  110 S(KLM)=S(KLM)+BIL*PABL
      GO TO 6
    4 IF (NVO.EQ.0) GO TO 5
C***  CALCULATE D2DJ/DXDX
      DO 332 KA=1,KK
      DO 332 L=KA,KK
      L1=IMG(KA)
      L2=IMG(L)
      PABL=ABL1(L1)*ABL1(L2)
      SUM=0.0
      DO 333 I=1,3
      SUM1=0.0
      DO 334 K=1,3
  334 SUM1=SUM1+G(I,K)*ALIN(K,KA)
  333 SUM=SUM+SUM1*ALIN(I,L)
      ABL2=(SUM-PABL)/D(J)
      IF (L1.GT.L2) GO TO 335
      KLM=INDEX(L1,L2)
      GO TO 336
  335 KLM=INDEX(L2,L1)
  336 S(KLM)=S(KLM)+BIL*(PABL-DELTAD(J)*ABL2)
  332 CONTINUE
    5 IF (IGIT.NE.1) GO TO 6
      DO 115 L1=1,NVG
C***  CALCULATE D2DJ/DADA
      LH1=IFELD(L1)
      LL=NVO+L1
      DO 112 L2=LL,NV1
      L3=L2-NVO
      LH2=IFELD(L3)
      SUM=0.0
      DO 111 K=1,3
      SUM=SUM+GABL2A(K,LH1,LH2)*PDELX(1,K)
      SUM=SUM+GABL2B(K,LH1,LH2)*PDELX(2,K)
  111 SUM=SUM+GABL2C(K,LH1,LH2)*PDELX(3,K)
      SUM=0.5*SUM
      PABL=ABL1(LL)*ABL1(L2)
      ABL2=(SUM-PABL)/D(J)
      IF (LL.GT.L2) GO TO 42
      KLM=INDEX(LL,L2)
      GO TO 43
   42 KLM=INDEX(L2,LL)
   43 S(KLM)=S(KLM)+BIL*(PABL-DELTAD(J)*ABL2)
  112 CONTINUE
      IF (NVO.EQ.0) GO TO 115
C***  CALCULATE D2DJ/DXDA
      DO 337 K=1,KK
      L2=IMG(K)
      PABL=ABL1(LL)*ABL1(L2)
      SUM=0.0
      DO 338 I=1,3
  338 SUM=SUM+GIKDAL(I,LH1)*ALIN(I,K)
      IF (L2.GT.LL) GO TO 339
      KLM=INDEX(L2,LL)
      GO TO 340
  339 KLM=INDEX(LL,L2)
  340 S(KLM)=S(KLM)+BIL*(PABL-DELTAD(J)*ABL2)
  337 CONTINUE
  115 CONTINUE
    6 IF (IDOB.EQ.0) GO TO 200
C***  CALCULATE MATRIX AND VECTOR FOR DISTANCE VARIABLES
      IG=LJ(J)
      HIL=RAT(IG)*RAT(IG)
      BIL=2.0*W(J)
      DO 118 I=1,NVD
      IF (LDR(IG).NE.I) GO TO 118
      K1=NV1+I
      T(K1)=T(K1)-BIL*DELTAD(J)*RAT(IG)
      DO 116 K=I,NVD
      IF (LDR(IG).NE.K) GO TO 116
      K2=NV1+K
      IF (K1.GT.K2) GO TO 46
      KLM=INDEX(K1,K2)
      GO TO 47
   46 KLM=INDEX(K2,K1)
   47 S(KLM)=S(KLM)+HIL*BIL
  116 CONTINUE
      DO 117 K=1,NV1
      IF (K.GT.K1) GO TO 48
      KLM=INDEX(K,K1)
      GO TO 117
   48 KLM=INDEX(K1,K)
  117 S(KLM)=S(KLM)-BIL*RAT(IG)*ABL1(K)
  118 CONTINUE
  200 CONTINUE
C***  CALCULATE R FACTORS
      HN=SQRT(RO)
      RWERT=HN/SQRT(RN)
      RSIGMA=HN/SQRT(FLOAT(M-NV))
      IF (NRE.EQ.0) GO TO 38
C***  MATRIX AND VECTOR IF LINEAR RESTRICTIONS ARE USED
      DO 133 L=1,NRE
      IF (LINH(L).EQ.1) GO TO 133
      L1=IZ1(L)
      DO 132 I=1,L1
      K1=KARI(L,I)
      PABL=2.0*WR(L)*RESKO(L,I)
      T(K1)=T(K1)-PABL*BEDING(L)
      DO 131 K=I,L1
      K2=KARI(L,K)
      IF (K1.LE.K2) GO TO 37
      L2=K1
      K1=K2
      K2=L2
   37 KLM=INDEX(K1,K2)
      S(KLM)=S(KLM)+PABL*RESKO(L,K)
  131 CONTINUE
  132 CONTINUE
  133 CONTINUE
   38 CONTINUE
C***  IS MATRIX SINGULAR...
      DO 123 I=1,NV
      K1=INDEX(I,I)
      IF (ABS(S(K1)).GT.1.0E-08) GO TO 123
      T(I)=0.0
      S(K1)=1.0
      I1=I+1
      IF (I1.GT.NV) GO TO 13
      DO 121 K=I1,NV
      K1=INDEX(I,K)
      S(K1)=0.0
  121 CONTINUE
   13 I1=I-1
      IF (I1.EQ.0) GO TO 14
      DO 122 K=1,I1
      K1=INDEX(K,I)
      S(K1)=0.0
  122 CONTINUE
   14 WRITE (NTOUT,1004) I
 1004 FORMAT (//11H***VARIABLE,I4,1X,63HCAUSES THE MATRIX TO BE SINGUL
     1AR OR DEGENERATE.FORCED INVARIANT)
  123 CONTINUE
      RETURN
      END
      SUBROUTINE MATOUT
C      PURPOSE... TO PUT OUT THE NORMAL OR CORRELATION MATRIX
      DIMENSION RCO(10), IX(3), NGIT(6), CHAR(2)
      REAL*8 ID,IX,KOR,CHAR,KOR1,KOR2,NAME,NGIT,NAME1,NAME2
C BK      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG                   0
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK4/ RN,RSIGMA,RWERT,ICVTST
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK9/ DELV(300),IFELD(6),IV(300),IV1(300),KVAR(3,200)
      COMMON S(50000)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
C BK      COMMON /BLOC18/ ROV,NC,NA,NS,IVF
      COMMON /BLOC18/ ROV,NC,NA,NS,IVF,IRNG1
      DATA IROW/'ROW '/
      DATA IX/' X',' Y',' Z'/
      DATA NGIT/' A',' B',' C','ALPHA','BETA','GAMMA'/
      DATA CHAR/'CELL',' '/
      INDEX(I,K)=NV*I+K-(I*(I-1))/2
      IF (MEIG.EQ.0) GO TO 4
      IZEIL=IZEIL+1
      CALL KOPF1
      WRITE (NTOUT,1006)
 1006 FORMAT ( 48H APPROXIMATE EIGENVALUES OF MATRIX (LAST CYCLE):/)
      GO TO 5
    4 LINE=NV/10+7
      IF (MATCOR.EQ.1) GO TO 1
      IF (IZEIL+LINE.GT.LETZTE) CALL KOPF
      IZEIL=IZEIL+5
      WRITE (NTOUT,1000)
 1000 FORMAT (//' NORMAL MATRIX'/)
      GO TO 2
    1 CALL RHO
      SGM=2*ROV/FLOAT(M+NRE-NRH-NV)
      SGM=SQRT(SGM)
      IF (IZEIL+10.GT.LETZTE) CALL KOPF
      IZEIL=IZEIL+4
      WRITE (NTOUT,2000)
 2000 FORMAT (/// 30H ESTIMATED STANDARD DEVIATIONS)
      IF (NVO.EQ.0) GO TO 202
      WRITE (NTOUT,2004)
 2004 FORMAT (/  5H ATOM,8X,  1HX,11X,  1HY,11X,  1HZ/)
      IZEIL=IZEIL+3
      DO 200 I=1,N1
      DO 201 K=1,3
      RCO(K)=0.0
      KV=KVAR(K,I)
      IF (KV.EQ.0) GO TO 201
      KLI=INDEX(KV,KV)
      RCO(K)=SQRT(S(KLI))*SGM
  201 CONTINUE
      CALL KOPF1
  200 WRITE (NTOUT,2001) ID(I),(RCO(J),J=1,3)
      IZEIL=IZEIL+3
 2001 FORMAT (1X,A6,1X,3(E10.4,2X))
  202 IF (NVG.EQ.0) GO TO 204
      WRITE (NTOUT,2005)
 2005 FORMAT (/7X,  1HA,11X,  1HB,11X,  1HC,8X,  5HALPHA,7X,  4HBETA,8X,
     1  5HGAMMA/)
      IZEIL=IZEIL+3
      SGM=1.4142*RSIGMA
      DO 205 I=1,6
  205 RCO(I)=0.0
      DO 206 L=1,NVG
      I=IFELD(L)
      L1=NVO+L
      KLI=INDEX(L1,L1)
  206 RCO(I)=SQRT(S(KLI))*SGM
      WRITE (NTOUT,2006) (RCO(I),I=1,6)
 2006 FORMAT (1X,6(E10.4,2X))
  204 IF (ICOR) 10,11,11
   10 IF (IZEIL+LINE.GT.LETZTE) CALL KOPF
      IZEIL=IZEIL+6
      WRITE (NTOUT,1001)
      GO TO 2
   11 IF (IZEIL+10.GT.LETZTE) CALL KOPF
      WRITE (NTOUT,1010) CORR
 1010 FORMAT (//' ELEMENTS OF CORRELATION MATRIX IN MAGNITUDE', 13H GRE
     1ATER THAN,F5.2/)
 1001 FORMAT (//' CORRELATION MATRIX'//)
      INONE=0
    2 DO 102 I=1,NV
      IF (ICOR.GT.0.AND.MATCOR.EQ.1) GO TO 12
      CALL KOPF1
      WRITE (NTOUT,1002) IROW,I
 1002 FORMAT (1X,A4,I4)
      K1=0
   12 DO 100 K=I,NV
      KLM=INDEX(I,K)
      IF (MATCOR.NE.1) GO TO 3
      KLI=INDEX(I,I)
      KLK=INDEX(K,K)
C***  CALCULATE CORRELATION MATRIX
      S(KLM)=S(KLM)/SQRT(S(KLI)*S(KLK))
      IF (ICOR.LT.0) GO TO 3
      IF (ABS(S(KLM)).LE.CORR.OR.I.EQ.K) GO TO 100
      INONE=1
      CALL KOPF1
      IL=I
  101 IF (IL.GT.NVO) GO TO 104
      NAME=ID(IV1(IL))
      KOR=IX(IV(IL))
      GO TO 107
  104 IFD=IFELD(IL-NVO)
      IF (IFD.GT.3) GO TO 105
      NAME=CHAR(1)
      KOR=NGIT(IFD)
      GO TO 107
  105 NAME=NGIT(IFD)
      KOR=CHAR(2)
  107 IF (IL.EQ.K) GO TO 108
      NAME1=NAME
      KOR1=KOR
      IL=K
      GO TO 101
  108 NAME2=NAME
      KOR2=KOR
      WRITE (NTOUT,1011) NAME1,KOR1,NAME2,KOR2,S(KLM)
 1011 FORMAT (1X,A6,A2,3X,  1H-,3X,A6,A2,8X,E10.4)
      GO TO 100
    3 K1=K1+1
      RCO(K1)=S(KLM)
      IF (K1.NE.10) GO TO 100
      K1=0
      CALL KOPF1
C***  PUT OUT MATRIX
      WRITE (NTOUT,1003) (RCO(J),J=1,10)
 1003 FORMAT (1X,10E11.4)
  100 CONTINUE
      IF (ICOR.GT.0.AND.MATCOR.EQ.1) GO TO 102
      IF (K1.EQ.0) GO TO 102
      CALL KOPF1
      WRITE (NTOUT,1003) (RCO(J),J=1,K1)
  102 CONTINUE
      IF (MATCOR.EQ.0) GO TO 203
      IF (ICOR.LE.0.OR.INONE.EQ.1) GO TO 203
      CALL KOPF1
      WRITE (NTOUT,2003)
 2003 FORMAT (  5H NONE)
  203 CALL KOPF1
      WRITE (NTOUT,1004)
 1004 FORMAT (1H )
      IF (MATCOR.NE.0) RETURN
      IZEIL=IZEIL+1
      CALL KOPF1
      WRITE (NTOUT,1005)
 1005 FORMAT ( 14H NORMAL VECTOR/)
    5 K1=0
C***  PUT OUT VECTOR OR APPR. EIGENVALUES OF NORMAL MATRIX
      DO 103 I=1,NV
      K1=K1+1
      IF (MEIG.EQ.0) GO TO 6
      RCO(K1)=DGEL(I)
      GO TO 7
    6 RCO(K1)=T(I)
    7 IF (K1.NE.10) GO TO 103
      K1=0
      CALL KOPF1
      WRITE (NTOUT,1003) (RCO(J),J=1,10)
  103 CONTINUE
      IF (K1.EQ.0) GO TO 8
      CALL KOPF1
      WRITE (NTOUT,1003) (RCO(J),J=1,K1)
      CALL KOPF1
      WRITE (NTOUT,1004)
    8 MEIG=0
      RETURN
      END
      SUBROUTINE INVERT
C      MATRIX INVERSION ROUTINE
C      -----------------------
C       BASED ON THE ALGORITHM OF BUSING AND LEVY,
C       COMM.ACM,5,445(1962)
C
C       INVERT IS CALLED BY...MAIN PROGRAM
C       AND CALLS...INCH.
C
C       IF NDIAG=1, A DIAGONAL MATRIX
C       APPROXIMATION IS RUN.
C
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON S(50000)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
      COMMON /BLOC11/ II
      INDEX(J1,J2)=J1*NV+J2-(J1*(J1-1))/2
      IF (NV.GT.1) GO TO 80
      S(1)=1.0/S(1)
      RETURN
   80 CONTINUE
      IF (NDIAG.NE.1) GO TO 2
      DO 99 I=1,NV
      K1=INDEX(I,I)
      S(K1)=1.0/S(K1)
      I1=I+1
      IF (I1.GT.NV) GO TO 99
      DO 98 K=I1,NV
      K1=INDEX(I,K)
      S(K1)=0.0
   98 CONTINUE
   99 CONTINUE
      RETURN
    2 DO 102 II=1,NV
      K1=INDEX(II-1,II)
      S(K1)=FLOAT(II)
      BIGS=0.0
      DO 100 J=II,NV
      K2=INDEX(J,J)
      PABL=ABS(S(K2))
      IF (PABL.LE.BIGS) GO TO 100
      BIGS=PABL
      S(K1)=FLOAT(J)
  100 CONTINUE
      CALL INCH
      K3=INDEX(II,II)
      PABL=ABS(S(K3))
      IF (PABL.GT.0.0) GO TO 1
      WRITE (NTOUT,1000) II
 1000 FORMAT (1HO//5H *** , 25HWARNING: DIAGONAL ELEMENT,I4,  9H VANISHE
     1S,4H ***//1H , 20HJOB KILLED BY INVERT)
      STOP
C***   STORE APPR. EIGENVALUES OF MATRIX S INTO VECTOR DGEL
    1 DGEL(II)=S(K3)
      S(K3)=1.0/S(K3)
      I=II+1
      IF (I.GT.NV) GO TO 10
      DO 101 K=I,NV
      K4=INDEX(II-1,K)
      K5=INDEX(II,K)
  101 S(K4)=-S(K3)*S(K5)
      DO 102 J=I,NV
      DO 102 K=J,NV
      K1=INDEX(J,K)
      K2=INDEX(II-1,K)
      K3=INDEX(II,J)
      S(K1)=S(K1)+S(K2)*S(K3)
  102 CONTINUE
   10 CONTINUE
      I1=NV-1
      LL=NV
      DO 106 L=1,I1
      II=LL-L
      I=II+1
      IF (I.GT.NV) GO TO 11
      DO 103 J=I,NV
      K1=INDEX(II,J)
  103 S(K1)=0.0
      DO 105 J=I,NV
      DO 104 K=J,NV
      K1=INDEX(II,K)
      K2=INDEX(II-1,J)
      K3=INDEX(J,K)
      S(K1)=S(K1)+S(K2)*S(K3)
      IF (J.EQ.K) GO TO 104
      K1=INDEX(II,J)
      K2=INDEX(II-1,K)
      S(K1)=S(K1)+S(K2)*S(K3)
  104 CONTINUE
      K1=INDEX(II,II)
      K2=INDEX(II-1,J)
      K3=INDEX(II,J)
      S(K1)=S(K1)+S(K2)*S(K3)
  105 CONTINUE
   11 CONTINUE
  106 CALL INCH
      RETURN
      END
      SUBROUTINE INCH
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON S(50000)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
      COMMON /BLOC11/ II
      INDEX(J1,J2)=J1*NV+J2-(J1*(J1-1))/2
      K1=INDEX(II-1,II)
      HN=S(K1)+0.5
      NH=INT(HN)
      IF (NH.EQ.II) RETURN
      DO 100 J=II,NV
      K1=INDEX(II,J)
      IF (J.LT.NH) GO TO 1
      K2=INDEX(NH,J)
      PABL=S(K1)
      S(K1)=S(K2)
      S(K2)=PABL
      GO TO 100
    1 K2=INDEX(J,NH)
      PABL=S(K1)
      S(K1)=S(K2)
      S(K2)=PABL
  100 CONTINUE
      K1=INDEX(II,II)
      K2=INDEX(II,NH)
      PABL=S(K1)
      S(K1)=S(K2)
      S(K2)=PABL
      RETURN
      END
      SUBROUTINE NEWPAR
C      PURPOSE... PARAMETERCHANGES= MATRIX * VECTOR.CORRECTS PARAMETERS.
      DIMENSION NN1(3), NN2(3), DELVG(300)
      REAL*8 ID
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK7/ IDEC(100),KA(2,150),LCELL(6)
      COMMON /BLOCK8/ SI(3,150),SIGN(3,3,150),SYMOPS(10,200),B(3,3,200),
     1BI(3,200),IND,MTR(200)
      COMMON /BLOCK9/ DELV(300),IFELD(6),IV(300),IV1(300),KVAR(3,200)
      COMMON S(50000)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
C BK      COMMON /BLOC13/ ISYS,CVGTST
      COMMON /BLOC13/ ISYS,CVGTST,LSYS
      COMMON /BLOC14/ IZ1(35),WR(35),BEDING(35),RESKO(35,20),KARI(35,20)
     1,LINH(35),LIND(35)
C BK      COMMON /BLOC18/ ROV,NC,NA,NS,IVF
      COMMON /BLOC18/ ROV,NC,NA,NS,IVF,IRNG1
      INDEX(I,K)=NV*I+K-(I*(I-1))/2
      NVF=0
      IBED=0
C***  CALCULATE PARAMETER CHANGES
      DO 103 I=1,NV
      SUM=0.0
      DO 101 K=I,NV
      K2=INDEX(I,K)
  101 SUM=SUM+S(K2)*T(K)
      K1=I-1
      IF (K1.EQ.0) GO TO 103
      DO 102 K=1,K1
      K2=INDEX(K,I)
  102 SUM=SUM+S(K2)*T(K)
  103 DELV(I)=SUM
      NN1(1)=NVO
      NN1(2)=NVG
      NN1(3)=NVD
      NN2(1)=0
      NN2(2)=NVO
      NN2(3)=NV1
      DD=0.001
C***  CHECK FUDGE FACTORS
      DO 105 L=1,3
      IF (FUDGE(L).LT.DD) GO TO 105
      K1=NN1(L)
      IF (K1.EQ.0) GO TO 105
      DO 104 I=1,K1
      K2=NN2(L)+I
  104 DELV(K2)=DELV(K2)*FUDGE(L)
  105 CONTINUE
      FK=1.0
      IF (IVF.LT.1) GO TO 108
      CALL RHO
      ROI=ROV
C***  STORE OLD PARAMETERS
  108 IF (NVO.EQ.0) GO TO 118
      DO 111 I=1,N1
      DO 111 K=1,3
  111 XOLD(K,I)=X(K,I)
  118 IF (IGIT.NE.1) GO TO 119
      DO 112 I=1,6
  112 AOLD(I)=A(I)
  119 IF (IDOB.NE.1) GO TO 400
      DO 110 I=1,NDO
  110 DOBOLD(I)=DOB(I)
  400 IF (NVO.EQ.0) GO TO 2
C***  NEW COORDINATES
      DO 106 I=1,NVO
      K1=IV1(I)
      K2=IV(I)
  106 X(K2,K1)=XOLD(K2,K1)+DELV(I)*FK
      IF (NRE.EQ.0) GO TO 116
      GO TO 113
      ENTRY LINPAR
      IBED=1
C***  CALCULATE TOTAL PARAMETER CHANGES: X(ACTUAL CYCLE)-X(START)
  113 DO 114 I=1,NVO
      K1=IV1(I)
      K2=IV(I)
  114 DELVG(I)=X(K2,K1)-XINI(K2,K1)
C***  CALCULATE LINEAR RESTRICTIONS OF CURRENT CYCLE
      DO 215 L=1,NRE
      L1=IZ1(L)
      BEDING(L)=0.
      DO 115 I=1,L1
      IF (I.EQ.LIND(L)) GO TO 115
      K1=KARI(L,I)
      BEDING(L)=BEDING(L)+RESKO(L,I)*DELVG(K1)
  115 CONTINUE
C***  COORDINATES DEPENDENT BY HARD RESTRICTIONS
      IF (LINH(L).EQ.0) GO TO 215
      L1=LIND(L)
      K2=MOD(KARI(L,L1),3)
      K1=1+(KARI(L,L1)-K2)/3
      X(K2,K1)=-BEDING(L)/RESKO(L,L1)+XINI(K2,K1)
      BEDING(L)=0.0
  215 CONTINUE
  116 DO 200 NL=1,N1
      IF (MTR(NL).EQ.0) GO TO 200
C***  IF ATOM IS IN SPECIAL POSITION
      DO 202 K=1,3
      XK=0.0
      DO 201 L=1,3
  201 XK=XK+B(K,L,NL)*X(L,NL)
  202 X(K,NL)=XK+BI(K,NL)
  200 CONTINUE
      IF (NEQU.EQ.0) GO TO 2
      DO 107 I=1,NEQU
      K1=KA(1,I)
      K2=KA(2,I)
      DO 107 K=1,3
      IF (IBED.EQ.1) GO TO 117
      XOLD(K,K2)=X(K,K2)
  117 X(K,K2)=SIGN(K,1,I)*X(1,K1)+SIGN(K,2,I)*X(2,K1)+SIGN(K,3,I)*X(3,K1
     1)+SI(K,I)
  107 CONTINUE
    2 IF (IBED.EQ.1) RETURN
      IF (IGIT.NE.1) GO TO 3
C***  NEW CELL PARAMETERS
      DO 109 I=1,NVG
      K=NVO+I
      K1=IFELD(I)
  109 A(K1)=AOLD(K1)+DELV(K)*FK
      GO TO (3,5,6,7), ISYS
C***   TRIGONAL,HEXAGONAL,TETRAGONAL
    5 A(2)=A(1)
      GO TO 3
C***   RHOMBOHEDRAL
    6 A(2)=A(1)
      A(3)=A(1)
      A(5)=A(4)
      A(6)=A(4)
      GO TO 3
C***   CUBIC
    7 A(2)=A(1)
      A(3)=A(1)
    3 IF (IDOB.NE.1) GO TO 121
C***  NEW PRESCRIBED DISTANCES
      DO 120 I=1,NDO
      IF (LDR(I).EQ.0) GO TO 120
      K=NV1+LDR(I)
      DOB(I)=DOBOLD(I)+DELV(K)*RAT(I)*FK
  120 CONTINUE
  121 NVF=NVF+1
      IF (NVF.GT.IVF) GO TO 122
      CALL RHO
      IF (ROV.LE.ROI) GO TO 122
      FK=FK*0.5
      GO TO 400
  122 IF (FK.GE.1.0) RETURN
      WRITE (NTOUT,1000) FK
 1000 FORMAT (' VARIABLE DAMPING FACTOR IS USED FOR COORDINATE CHANGES.'
     1,' THE FACTOR APPLIED WAS ',F7.4)
      IZEIL=IZEIL+1
      CALL KOPF1
      RETURN
      END
      SUBROUTINE DISOUT
C      PURPOSE... TO PUT OUT DISTANCES
      INTEGER BT
	REAL*8 ID
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOC15/ NAPID,IRNG,KRNG
      COMMON /BLOC16/ ITYPE(200),BT(15),IA,NZA,ICON(40,11),MD(40,18),NCA
     1PID,CAPID(15,4),IAT,IAB,IBT(40,4),DOT,WEIGHT(15,3),DELTAX(3),DELTA
     20(3),DA,DB,ANG,TOT(40,4),OTO(40,6),PI180,IMOV(3),MSYMAB(200),DOBIN
     3(400),KTYPE(20),IDKTYP
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
C BK     1IDIMM
      IZEIL=IZEIL+2
      WRITE (NTOUT,1002)
      IF (IZEIL+4.GT.LETZTE) CALL KOPF
      WRITE (NTOUT,1000) NCC
 1000 FORMAT (35H INTERATOMIC DISTANCES BEFORE CYCLE,I4)
      IF (NCAPID.NE.0) WRITE (NTOUT,1005) NCAPID
 1005 FORMAT (  1H+,39X, 13HOF APID CYCLE,I3)
      IZEIL=IZEIL+4
      WRITE (NTOUT,1003)
 1003 FORMAT (  //' ATOM 1 ATOM 2 CALCULATED D PRESCRIBED DO',
	15X,'DO-D    WEIGHT     W*(DO-D)')
      DO 100 J=1,M
      CALL KOPF1
      K1=NN(J,1)
      K2=NN(J,2)
      K=LJ(J)
      SQRTW=SQRT(W(J))
      DW=SQRTW*DELTAD(J)
      WRITE (NTOUT,1001) ID(K1),ID(K2),D(J),DOB(K),DELTAD(J),SQRTW,DW
 1001 FORMAT (2(1X,A6),F10.4,F14.4,F12.4,F10.4,F12.4)
  100 CONTINUE
      WRITE (NTOUT,1002)
 1002 FORMAT (//)
      IZEIL=IZEIL+1
      CALL KOPF1
      RETURN
      END
      SUBROUTINE KOPF
C     PURPOSE... COUNTS PAGES, WRITES TITLE.
C     REAL*8 ID,JNAME,JDATE
C     INTEGER*2 DATE(3)
      character*9 dte
      REAL*8 ID
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
C     EQUIVALENCE (JDATE,DATE(1))
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
C     CALL JNMDAT(JNAME,JDATE)
      call xrsdate(DTE)
      IPAGE=IPAGE+1
      IZEIL=0
      LETZTE=55
C
C REFORMATED TO AVOID HORRIBLE OUTPUT IN NOTEPAD ON PC'S RJA VTX JAN02
      WRITE (NTOUT,1000) ITIT,DTE,IPAGE
 1000 FORMAT (///'   DLS-76-02',7X,18A4,'  DATE: ',A9,7X,'PAGE',I4
     1        ////)
      RETURN
      END
      SUBROUTINE KOPF1
C      PURPOSE... LINE COUNTER.
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      IZEIL=IZEIL+1
      IF (IZEIL.LE.LETZTE) RETURN
      CALL KOPF
      RETURN
      END
      SUBROUTINE PAROUT
C      PURPOSE...  TO PUT OUT THE OLD AND NEW PARAMETERS
      DIMENSION KXXK(3), A1(6), A1OLD(6), ITEXT(2,2), DXG(3)
      DIMENSION XXN(2), XXD(2), XXO(2), XXDA(2), DSH(2)
      REAL*8 ID,ITEXT
      REAL IHR,IFAK,LCELL
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK4/ RN,RSIGMA,RWERT,ICVTST
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK7/ IDEC(100),KA(2,150),LCELL(6)
      COMMON /BLOCK9/ DELV(300),IFELD(6),IV(300),IV1(300),KVAR(3,200)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
      COMMON /BLOC14/ IZ1(35),WR(35),BEDING(35),RESKO(35,20),KARI(35,20)
     1,LINH(35),LIND(35)
      COMMON /BLOC15/ NAPID,IRNG,KRNG
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
C BK     1IDIMM
      COMMON /BLOC18/ ROV,NC,NA,NS,IVF,IRNG1
      DATA KXXK/'X   ','Y   ','Z   '/,KATOM/'ATOM'/
      DATA ITEXT/'LAST CYC','THIS TRA','LE      ','NSLATION'/
      IF (ICVTST.LT.2) GO TO 200
      IF (IGIT.EQ.1) GO TO 19
      DO 20 I=1,N1
      DO 20 K=1,3
      IF (IFAK(K,I).NE.IHR) GO TO 20
      IF (ABS(XINI(K,I)).GT.1.E-10) GO TO 21
   20 CONTINUE
   19 IF (ICVTST.EQ.3) RETURN
      CALL PAREX
      RETURN
   21 CALL KOPF
      WRITE (NTOUT,1021) ITEXT(ICVTST-1,1),ITEXT(ICVTST-1,2)
 1021 FORMAT ( 30H TOTAL PARAMETER SHIFTS AFTER ,2A8 //2(4X, 10H PARAMET
     1ER,4X,  7HINITIAL,5X,  6HCHANGE,6X,  5HFINAL,6X,  5HSHIFT,3X)/)
      IZEIL=IZEIL+5
      GO TO 2
C*** OUTPUT IF RANDOM COORDINATES ARE USED
  200 IF (IRNG.EQ.0) GO TO 7
      IF (NCYC.EQ.0) GO TO 30
      IF (NCYC.EQ.1.AND.IZEIL.GT.0) CALL KOPF
      IZEIL=IZEIL+3
      CALL RHO
      WRITE (NTOUT,1013) KRNG,NCYC,RWERT,RSIGMA
 1013 FORMAT (//' SET',I4,1H:,2X, 18HAFTER CYCLE NUMBER,I5,11X,  3HR =,F
     113.8,3X,6HSIGMA=,F13.8)
      IF (NRE.EQ.NRH) GO TO 96
      WRITE (NTOUT,1026) ROV
 1026 FORMAT (  1H+,85X,  4HRHO=,F13.8)
      GO TO 96
   30 IF (IZEIL.GT.0) CALL KOPF
      WRITE (NTOUT,1030) KRNG,IRNG1
 1030 FORMAT (//' RANDOM COORDINATES (STARTING VALUES) OF SET',I3/ 47H -
     1---------------------------------------------//30H GENERATOR START
     2ING NUMBER WAS,I12)
      IZEIL=IZEIL+4
   96 WRITE (NTOUT,1027)
 1027 FORMAT (  1H )
      DO 97 I=1,N
      CALL KOPF1
   97 WRITE (NTOUT,1014) KATOM,ID(I),(X(K,I),K=1,3)
 1014 FORMAT (1X,A4,1X,A6,3F10.5)
      RETURN
C*** END OUTPUT OF RANDOM COORDINATES
    7 IF(IGIT.EQ.0) GO TO 8
      PI180=180./3.141592654
      DO 98 K=1,3
      A1(K)=A(K)
   98 A1OLD(K)=AOLD(K)
      DO 99 K=4,6
      A1(K)=A(K)*PI180
   99 A1OLD(K)=AOLD(K)*PI180
    8 IZEIL=IZEIL+8
      WRITE (NTOUT,1025)
 1025 FORMAT (///)
      IF (IZEIL+4.LT.LETZTE) GO TO 90
      CALL KOPF
      IZEIL=IZEIL+5
   90 WRITE (NTOUT,1001) NCYC
 1001 FORMAT (//' PARAMETERS AFTER CYCLE ',I4)
      WRITE (NTOUT,1002)
 1002 FORMAT (/2(4X, 10H PARAMETER,6X,  3HOLD,7X,  6HCHANGE,6X,  3HNEW,1
     16X)/)
      IF (IGIT.EQ.0) GO TO 2
      K=NVO
C*** PUT OUT CELL PARAMETERS
      DO 101 I=1,6
      CALL KOPF1
      XX1=A1(I)-A1OLD(I)
  101 WRITE (NTOUT,1003) LCELL(I),A1OLD(I),XX1,A1(I)
 1003 FORMAT (2(8X,A4,2X,3(1X,F10.5),14X))
      WRITE (NTOUT,1003)
      CALL KOPF1
    2 IF (NVO.EQ.0) GO TO 3
C*** PUT OUT ATOM PARAMETERS
      I98=-1
   10 I98=I98+2
      I89=I98
      I99=I98+1
      IF (IZEIL+5.GT.LETZTE) CALL KOPF
      CALL KOPF1
      IF (I99.LE.N1) GO TO 12
      I90=I98
      I91=1
      GO TO 13
   12 I90=I99
      I91=2
   13 WRITE (NTOUT,1004) (KATOM,ID(I),I=I89,I90)
 1004 FORMAT (1X,A4,1X,A6,50X,A4,1X,A6)
      IF (ICVTST.LT.2) GO TO 201
      DO 202 K=1,3
      CALL KOPF1
      XXN(1)=X(K,I89)
      XXO(1)=XINI(K,I89)
      XXD(1)=XXN(1)-XXO(1)
      XXDA(1)=XXD(1)*A(K)
      IF (I91.EQ.1) GO TO 202
      XXN(2)=X(K,I90)
      XXO(2)=XINI(K,I90)
      XXD(2)=XXN(2)-XXO(2)
      XXDA(2)=XXD(2)*A(K)
  202 WRITE (NTOUT,1022) (KXXK(K),XXO(I),XXD(I),XXN(I),XXDA(I),I=1,I91)
 1022 FORMAT (2(8X,A4,2X,4(1X,F10.5),3X))
      DSH1=0.0
      DSH2=0.0
      DO 203 I=1,3
      DO 203 K=1,3
      DSH1=DSH1+G(I,K)*(X(I,I89)-XINI(I,I89))*(X(K,I89)-XINI(K,I89))
      IF (I91.EQ.1) GO TO 203
      DSH2=DSH2+G(I,K)*(X(I,I90)-XINI(I,I90))*(X(K,I90)-XINI(K,I90))
  203 CONTINUE
      DSH(1)=SQRT(DSH1)
      IF (I91.EQ.2) DSH(2)=SQRT(DSH2)
      WRITE (NTOUT,1023) (DSH(I),I=1,I91)
 1023 FORMAT (2(8X,'MAGNITUDE',31X,F10.5,3X))
      CALL KOPF1
      GO TO 204
  201 DO 102 K=1,3
      CALL KOPF1
      XXN(1)=X(K,I89)
      XXO(1)=XOLD(K,I89)
      XXD(1)=XXN(1)-XXO(1)
      IF (I91.EQ.1) GO TO 102
      XXN(2)=X(K,I90)
      XXO(2)=XOLD(K,I90)
      XXD(2)=XXN(2)-XXO(2)
  102 WRITE (NTOUT,1003) (KXXK(K),XXO(I),XXD(I),XXN(I),I=1,I91)
  204 IF (I99.LT.N1) GO TO 10
      IF (ICVTST.LT.2) GO TO 210
      IF (IRNG.GT.0) GO TO 3
      RETURN
  210 IF (NEQU.EQ.0) GO TO 3
      IF (IZEIL+6+NEQU.GT.LETZTE) CALL KOPF
      IZEIL=IZEIL+6
      WRITE (NTOUT,1015)
 1015 FORMAT (//' SYMMETRICALLY DEPENDENT ATOMS'/)
      WRITE (NTOUT,1016)
 1016 FORMAT (1X,4HATOM,6X,1HX,8X,1HY,8X,1HZ/)
      NH=N1+1
      DO 106 I=NH,N
      CALL KOPF1
  106 WRITE (NTOUT,1017) ID(I),(X(K,I),K=1,3)
 1017 FORMAT (1X,A6,3(F8.5,1X))
    3 IF (NRE.EQ.0) GO TO 6
C*** PUT OUT LINEAR RESTRICTIONS
      IF (IZEIL+6+N1.GT.LETZTE) CALL KOPF
      IZEIL=IZEIL+6
      WRITE (NTOUT,1007)
 1007 FORMAT (//' DIFFERENCE VECTORS TO INITIAL COORDINATES'/)
      WRITE (NTOUT,1008)
 1008 FORMAT (1X,4HATOM,11X,2HDX,9X,2HDY,9X,2HDZ/)
      DO 104 I=1,N1
      CALL KOPF1
      DO 108 K=1,3
  108 DXG(K)=X(K,I)-XINI(K,I)
  104 WRITE (NTOUT,1009) ID(I),(DXG(K),K=1,3)
 1009 FORMAT (1X,A6,2X,3F11.5)
      IF (NRE.EQ.NRH) GO TO 6
      IF (IZEIL+8+NRE-NRH.GE.LETZTE) CALL KOPF
      WRITE (NTOUT,1010) NCYC
 1010 FORMAT (//' LINEAR RESTRICTIONS AFTER CYCLE ',I4)
      IZEIL=IZEIL+6
      IF (NRH.EQ.0) GO TO 107
      WRITE (NTOUT,1018)
 1018 FORMAT (//' (HARD RESTRICTIONS ARE OMITTED FROM THIS LIST)')
      IZEIL=IZEIL+2
  107 WRITE (NTOUT,1011)
 1011 FORMAT (' NO OF CONDITION',9X,1HC,12X,8HWEIGHT*C,2X,20HWEIGHT*C/
     1NO OF TERMS/)
      DO 105 L=1,NRE
      IF (LINH(L).EQ.1) GO TO 105
      XX1=BEDING(L)*SQRT(WR(L))/WRF
      XX2=XX1/FLOAT(IZ1(L))
      CALL KOPF1
      WRITE (NTOUT,1012) L,BEDING(L),XX1,XX2
  105 CONTINUE
 1012 FORMAT (5X,I4,4X,3F16.5)
    6 IF (IDOB.EQ.0) RETURN
C***    PUT OUT DISTANCE VARIABLES
      IZEIL=IZEIL+3
      WRITE (NTOUT,1005)
 1005 FORMAT (//' PRESCRIBED DISTANCES'/  //7X,  3HOLD,7X,  6HCHANGE,
     16X,  3HNEW/)
      DO 103 I=1,NDO
      CALL KOPF1
      XX1=DOB(I)-DOBOLD(I)
  103 WRITE (NTOUT,1006) DOBOLD(I),XX1,DOB(I)
 1006 FORMAT (1X,3F11.5)
      RETURN
      END
      SUBROUTINE APID
C     PURPOSE... ADJUST PRESCRIBED INTERATOMIC DISTANCES ACCORDING TO
C     SUPPLIED FUNCTION
      INTEGER BT
      DIMENSION BLSTAR(2)
      REAL*8 ID
      REAL ISYMB
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOC12/ SYMSIG(3,3,40),KATOR(400),NN11(400,2)
C BK      COMMON /BLOC13/ ISYS,CVGTST
      COMMON /BLOC13/ ISYS,CVGTST,LSYS
      COMMON /BLOC16/ ITYPE(200),BT(15),IA,NZA,ICON(40,11),MD(40,18),NCA
     1PID,CAPID(15,4),IAT,IAB,IBT(40,4),DOT,WEIGHT(15,3),DELTAX(3),DELTA
     20(3),DA,DB,ANG,TOT(40,4),OTO(40,6),PI180,IMOV(3),MSYMAB(200),DOBIN
     3(400),KTYPE(20),IDK
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
C BK     1IDIMM
      COMMON /BLOC19/ RDD,ROTOL,ROTOU,RTOTL,RTOTU
      DATA BLSTAR(1),BLSTAR(2),CHAR/'  ',' *','C'/
      NCAPID=NCAPID+1
C
C     RESET WEIGHTS OF OUT OF RANGE DISTANCES (ANGLES)
C
C
C     CALCULATE NEW D(CALC)
C
      DO 4 J=1,M
      IF (MW(J).EQ.1) GO TO 3
      MW(J)=1
      W(J)=W(J)/WFAC
    3 DOBOLD(J)=DOB(J)
      IAT=NN(J,1)
      IAB=NN(J,2)
      CALL DIST
    4 D(J)=DOT
C     CALCULATE ANGLE AT BRIDGING ATOMS
      DO 35 NZ=1,NZA
      ITOT=ICON(NZ,10)
      I=0
      DO 35 I1=1,ITOT
   30 I=I+1
      IF (ICON(NZ,I+5).LE.0) GO TO 30
      DA=D(MD(NZ,I))
      DB=D(MD(NZ,I+10))
      DO 36 K=1,3
      DELTAX(K)=X(K,ICON(NZ,1))-X(K,ICON(NZ,I+1))
   36 DELTA0(K)=X(K,ICON(NZ,I+5))-X(K,ICON(NZ,I+1))
      CALL ANGLE
   35 TOT(NZ,I)=ANG
C     CALCULATE ANGLES AT CENTRAL ATOMS
      DO 40 NZ=1,NZA
      KK=0
      ITO=ICON(NZ,11)
      ITO1=ITO+1
      DO 40 I=2,ITO
      II=I+1
      DO 40 L=II,ITO1
      KK=KK+1
      DA=D(MD(NZ,I-1))
      DB=D(MD(NZ,L-1))
      DO 41 K=1,3
      DELTAX(K)=X(K,ICON(NZ,I))-X(K,ICON(NZ,1))
   41 DELTA0(K)=X(K,ICON(NZ,L))-X(K,ICON(NZ,1))
      CALL ANGLE
      OTO(NZ,KK)=ANG
      IF (ANG.GT.ROTOL.AND.ANG.LT.ROTOU) GO TO 40
      KLT=MD(NZ,KK+4)
      MW(KLT)=2
      W(KLT)=W(KLT)*WFAC
   40 CONTINUE
C     CALCULATE NEW PRESCRIBED VALUES FOR T-O DISTANCES
      DO 55 NZ=1,NZA
      ITOT=ICON(NZ,10)
      I=0
      DO 55 I1=1,ITOT
   45 I=I+1
      IF (ICON(NZ,I+5).LE.0) GO TO 45
      KL=MD(NZ,I)
      IF (ISYMB(KL).EQ.CHAR) GO TO 55
      II=IBT(NZ,I)
      TT=TOT(NZ,I)
      DTOT=TT-CAPID(II,4)
      IF (TT.GT.RTOTL.AND.TT.LT.RTOTU) GO TO 50
      KLT=MD(NZ,I+14)
      IF (MW(KLT).GT.1) GO TO 49
      MW(KLT)=2
      W(KLT)=W(KLT)*WFAC
   49 IF (ABS(WFAC-1.0).GT.1.E-10) GO TO 55
   50 DOB(KL)=CAPID(II,1)+CAPID(II,2)*DTOT+CAPID(II,3)*DTOT*DTOT
      DIDO=ABS(DOB(KL)-D(KL))/DOB(KL)
      IF (DIDO.LT.RDD) GO TO 55
      IF (MW(KL).GT.1) GO TO 55
      MW(KL)=2
      W(KL)=W(KL)*WFAC
   55 CONTINUE
C     CALCULATE NEW PRESCRIBED VALUES FOR O-O DISTANCES
      DO 60 NZ=1,NZA
      K=4
      ITO=ICON(NZ,11)
      ITO1=ITO+1
      DO 60 I=2,ITO
      II=I+1
      DO 60 L=II,ITO1
      K=K+1
      KL=MD(NZ,K)
      IF (ISYMB(KL).EQ.CHAR) GO TO 60
      DA=DOB(MD(NZ,I-1))
      DB=DOB(MD(NZ,L-1))
      DOB(KL)=SQRT(DA*DA+DB*DB+2./3.*DA*DB)
   60 CONTINUE
C     CALCULATE NEW PRESCRIBED VALUES FOR T-T DISTANCES
      DO 61 NZ=1,NZA
      ITOT=ICON(NZ,10)
      I=0
      DO 61 I1=1,ITOT
   62 I=I+1
      IF (ICON(NZ,I+5).LE.0) GO TO 62
      KL=MD(NZ,I)
      IF (ISYMB(KL).EQ.CHAR) GO TO 61
      DA=DOB(MD(NZ,I))
      DB=DOB(MD(NZ,I+10))
      CO=COS(CAPID(IBT(NZ,I),4)*PI180)
      DOB(MD(NZ,I+14))=SQRT(DA*DA+DB*DB-2.*DA*DB*CO)
   61 CONTINUE
      IF (IZEIL.GT.0) CALL KOPF
      IZEIL=IZEIL+4
      WRITE (NTOUT,1030) NCAPID
 1030 FORMAT (// 55H INTERATOMIC DISTANCES AND BOND ANGLES AFTER APID CY
     1CLE,I4/)
      CALL KOPF1
      WRITE (NTOUT,1032)
 1032 FORMAT (  7H ATOM 1,4X,  6HATOM 2,4X, 53HBOND  TYPE    OBSERVED D
     1   OLD DO   CHANGE    NEW DO,2X, 12HDO-DO(START),11X,  4HBOND,13X,
     2  5HANGLE)
      DO 65 NZ=1,NZA
      ITOT=ICON(NZ,10)
      IZEIL=IZEIL+1
      CALL KOPF1
      WRITE (NTOUT,1045)
      IF (IZEIL+4.GE.LETZTE) CALL KOPF
      L=0
      DO 67 L1=1,ITOT
   66 L=L+1
      IF (ICON(NZ,L+5).LE.0) GO TO 66
      KL=MD(NZ,L)
      CDD=DOB(KL)-DOBOLD(KL)
      DD=DOB(KL)-DOBIN(KL)
      NT1=ICON(NZ,1)
      NT2=ICON(NZ,L+1)
      NT3=ICON(NZ,L+5)
      KT2=MSYMAB(NT2)
      KT3=MSYMAB(NT3)
      IF (ISYMB(KL).EQ.CHAR) GO TO 63
      WRITE (NTOUT,1036) ID(NT1),ID(NT2),KTYPE(ITYPE(NT1)),KTYPE(ITYPE(K
     1T2)),KTYPE(ITYPE(KT3)),D(KL),BLSTAR(MW(KL)),DOBOLD(KL),CDD,DOB(KL)
     2,DD,ID(NT1),ID(NT2),ID(NT3),TOT(NZ,L)
 1036 FORMAT  (1X,2(A6,4X),3A4,F10.4,A2,F10.4,F9.4,1X,F9.4,F11.4,6X,A6,
     1         2H- ,A6,2H- ,A6,F8.1)
      GO TO 67
   63 WRITE (NTOUT,1033) ID(NT1),ID(NT2),KTYPE(ITYPE(NT1)),KTYPE(ITYPE(K
     1T2)),KTYPE(ITYPE(KT3)),D(KL),BLSTAR(MW(KL)),DOBOLD(KL),DOB(KL)
     2,DD,ID(NT1),ID(NT2),ID(NT3),TOT(NZ,L)
 1033 FORMAT (1X,2(A6,4X),3A4,F10.4,A2,F10.4,'  CONSTANT',F9.4,F11.4,6X,
     1         A6,2H- ,A6,2H- ,A6,F8.1)
   67 CALL KOPF1
      WRITE (NTOUT,1031)
      CALL KOPF1
      IF (IZEIL+6.GE.LETZTE) CALL KOPF
      K=4
      ITO=ICON(NZ,11)
      ITO1=ITO+1
      NT1=ICON(NZ,1)
      DO 68 I=2,ITO
      II=I+1
      DO 68 L=II,ITO1
      K=K+1
      KL=MD(NZ,K)
      CDD=DOB(KL)-DOBOLD(KL)
      DD=DOB(KL)-DOBIN(KL)
      NT2=ICON(NZ,I)
      NT3=ICON(NZ,L)
      K4=K-4
      IF (ISYMB(KL).EQ.CHAR) GO TO 64
      WRITE (NTOUT,1037) ID(NT2),ID(NT3),D(KL),BLSTAR(MW(KL)),DOBOLD(KL)
     1,CDD,DOB(KL),DD,ID(NT2),ID(NT1),ID(NT3),OTO(NZ,K4)
 1037 FORMAT (1X,2(A6,4X),11X,F11.4,A2,F10.4,F9.4,1X,F9.4,F11.4,6X,A6,
     1        2H- ,A6,2H- ,A6,F8.1)
      GO TO 74
   64 WRITE (NTOUT,1034) ID(NT2),ID(NT3),D(KL),BLSTAR(MW(KL)),DOBOLD(KL)
     1,DOB(KL),DD,ID(NT2),ID(NT1),ID(NT3),OTO(NZ,K4)
 1034 FORMAT (1X,2(A6,4X),11X,F11.4,A2,F10.4,'  CONSTANT',F9.4,F11.4,6X,
     1        A6,2H- ,A6,2H- ,A6,F8.1)
   74 CALL KOPF1
   68 CONTINUE
      WRITE (NTOUT,1031)
      CALL KOPF1
      IF (IZEIL+4.GE.LETZTE) CALL KOPF
      L=0
      DO 65 L1=1,ITOT
   69 L=L+1
      IF (ICON(NZ,L+5).LE.0) GO TO 69
      KL=MD(NZ,L+14)
      CDD=DOB(KL)-DOBOLD(KL)
      DD=DOB(KL)-DOBIN(KL)
      L5=L+5
      IF (ISYMB(KL).EQ.CHAR) GO TO 73
      WRITE (NTOUT,1038) ID(ICON(NZ,1)),ID(ICON(NZ,L5)),D(KL),BLSTAR(MW
     1(KL)),DOBOLD(KL),CDD,DOB(KL),DD
 1038 FORMAT (1X,2(A6,4X),11X,F11.4,A2,F10.4,F9.4,1X,F9.4,F11.4)
      GO TO 65
   73 WRITE (NTOUT,1035) ID(ICON(NZ,1)),ID(ICON(NZ,L5)),D(KL),BLSTAR(MW
     1(KL)),DOBOLD(KL),DOB(KL),DD
 1035 FORMAT (1X,2(A6,4X),11X,F11.4,A2,F10.4,'  CONSTANT',F9.4,F11.4)
   65 CALL KOPF1
 1031 FORMAT (1H )
      DO 71 J=1,M
      IF (MW(J).EQ.2) GO TO 72
   71 CONTINUE
      GO TO 70
   72 WRITE (NTOUT,1055) RDD,ROTOL,ROTOU,RTOTL,RTOTU
 1055 FORMAT (//' THE DISTANCES MARKED WITH AN ASTERISK ARE OUTSIDE', 22
     1H THE FOLLOWING RANGES:/ 30H FOR T-O BONDS: ABS(D-DO)/DO @,F5.2/ 2
     28H FOR O-O BONDS: O-T-O ANGLE:,F6.1,  2H -,F6.1/ 28H FOR T-T BONDS
     3: T-O-T ANGLE:,F6.1,  2H -,F6.1)
      IZEIL=IZEIL+4
      IF (WFAC.EQ.1.0) GO TO 70
      SQWF=SQRT(WFAC)
      WRITE (NTOUT,1060) SQWF
 1060 FORMAT ( 38H THEIR WEIGHTS HAVE BEEN MULTIPLIED BY,F10.1)
      IZEIL=IZEIL+1
   70 WRITE (NTOUT,1040)
 1040 FORMAT (///)
 1045 FORMAT (/)
      IZEIL=IZEIL+3
      CALL KOPF1
      RETURN
      END
      SUBROUTINE DIST
C     PURPOSE... CALCULATES DISTANCE BETWEEN TWO ATOMS
      DIMENSION DELTAX(3)
      INTEGER BT
      REAL*8 ID
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
      COMMON /BLOC16/ ITYPE(200),BT(15),IA,NZA,ICON(40,11),MD(40,18),NCA
     1PID,CAPID(15,4),IAT,IAB,IBT(40,4),DOT,WEIGHT(15,3),DELTA0(3),DELTA
     21(3),DA,DB,ANG,TOT(40,4),OTO(40,6),PI180,IMOV(3),MSYMAB(200),DOBIN
     3(400),KTYPE(20),IDKTYP
      DO 1 K=1,3
      DELTAX(K)=X(K,IAB)-X(K,IAT)
    1 CONTINUE
      DOTQ=0.0
      DO 2 J=1,3
      DO 2 K=1,3
      DOTQ=DOTQ+G(J,K)*DELTAX(J)*DELTAX(K)
    2 CONTINUE
      DOT=SQRT(DOTQ)
      RETURN
      END
      SUBROUTINE ANGLE
C     PURPOSE... CALCULATES ANGLES AT BRIDGING AND CENTRAL ATOMS
      INTEGER BT
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
      COMMON /BLOC16/ ITYPE(200),BT(15),IA,NZA,ICON(40,11),MD(40,18),NCA
     1PID,CAPID(15,4),IAT,IAB,IBT(40,4),DOT,WEIGHT(15,3),DELTAX(3),DELTA
     20(3),DA,DB,ANG,TOT(40,4),OTO(40,6),PI180,IMOV(3),MSYMAB(200),DOBIN
     3(400),KTYPE(20),IDKTYP
      SCAPR=0.
      DO 1 J=1,3
      DO 1 K=1,3
    1 SCAPR=SCAPR+G(J,K)*DELTAX(J)*DELTA0(K)
      COSA=SCAPR/DA/DB
      IF (ABS(COSA).GE.1.0) COSA=SIGN(1.0,COSA)
      ANG=ACOS(COSA)/PI180
      RETURN
      END
      SUBROUTINE EQUI
C     PURPOSE... TESTS IF DISTANCES NO JS AND NO MS ARE
C     SYMMETRICALLY EQUIVALENT
      DIMENSION NE(4), XA(3,2)
      REAL*8 ID
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK8/ SI(3,150),SIGN(3,3,150),SYMOPS(10,200),B(3,3,200),
     1BI(3,200),IND,MTR(200)
      COMMON /BLOC12/ SYMSIG(3,3,40),KATOR(400),NN11(400,2)
      COMMON /BLOC16/ ITYPE(200),BT(15),IA,NZA,ICON(40,11),MD(40,18),NCA
     1PID,CAPID(15,4),IAT,IAB,IBT(40,4),DOT,WEIGHT(15,3),DELTAX(3),DELTA
     20(3),DA,DB,ANG,TOT(40,4),OTO(40,6),PI180,IMOV(3),MSYMAB(200),DOBIN
     3(400),KTYPE(20),IDKTYP
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
C BK     1IDIMM
C     IF NOT THE SAME INDEPENDENT ATOMS
      IF ((NN11(JS,1).NE.NN11(MS,1).OR.NN11(JS,2).NE.NN11(MS,2)).AND.(NN
     111(JS,1).NE.NN11(MS,2).OR.NN11(JS,2).NE.NN11(MS,1))) GO TO 10
      DD=5.E-5
C     COMPARE DISTANCES
      IAT=NN(JS,1)
      IAB=NN(JS,2)
      CALL DIST
      DJS=DOT
      IAT=NN(MS,1)
      IAB=NN(MS,2)
      CALL DIST
      DMS=DOT
      IF (ABS(DMS-DJS).GT.DD) GO TO 10
C     PRODUCE ARTIFICIAL SHIFTS
      LL=2
      IF (NN11(JS,1).EQ.NN11(JS,2)) LL=1
      DO 1 L=1,LL
      DO 1 I=1,3
      XI=0.0
      DO 2 K=1,3
    2 XI=XI+B(I,K,NN11(JS,L))*SIN(FLOAT(K+3*L-3))
    1 XA(I,L)=XI+BI(I,NN11(JS,L))
      IF (LL.EQ.2) GO TO 3
      DO 4 I=1,3
    4 XA(I,2)=XA(I,1)
C     COORDINATES OF DEPENDENT ATOMS
    3 DO 5 LK=1,3,2
      DO 5 L=1,2
      IF (LK.EQ.1) NE(L)=NN(JS,L)-N1
      IF (LK.EQ.3) NE(L+2)=NN(MS,L)-N1
      IF (LK.EQ.3.AND.NN11(JS,1).NE.NN11(MS,1)) NE(L+2)=NN(MS,3-L)-N1
      IF (NE(L+LK-1).LE.0) GO TO 6
      DO 7 I=1,3
      XI=0.0
      DO 8 K=1,3
    8 XI=XI+SIGN(I,K,NE(L+LK-1))*XA(K,L)
    7 X(I,N+L+LK-1)=XI+SI(I,NE(L+LK-1))
      GO TO 5
    6 DO 9 I=1,3
    9 X(I,N+L+LK-1)=XA(I,L)
    5 CONTINUE
C     COMPARE DISTANCES AGAIN
      IAT=N+1
      IAB=N+2
      CALL DIST
      DJS=DOT
      IAT=N+3
      IAB=N+4
      CALL DIST
      DMS=DOT
      IF (ABS(DMS-DJS).GT.DD) GO TO 10
      ITEST=1
      RETURN
   10 ITEST=0
      RETURN
      END
      SUBROUTINE SETUP
C
C     PURPOSE: GENERATES ALL NONEQUIVALENT DISTANCES FROM SUPPLIED
C     CONNECTECTIVIY SPECIFICATIONS
C
      INTEGER BT
      REAL*8 ID
      REAL ISYMB
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK4/ RN,RSIGMA,RWERT,ICVTST
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOC12/ SYMSIG(3,3,40),KATOR(400),NN11(400,2)
C BK      COMMON /BLOC13/ ISYS,CVGTST
      COMMON /BLOC13/ ISYS,CVGTST,LSYS
      COMMON /BLOC16/ ITYPE(200),BT(15),IA,NZA,ICON(40,11),MD(40,18),NCA
     1PID,CAPID(15,4),IAT,IAB,IBT(40,4),DOT,WEIGHT(15,3),DELTAX(3),DELTA
     20(3),DA,DB,ANG,TOT(40,4),OTO(40,6),PI180,IMOV(3),MSYMAB(200),DOBIN
     3(400),KTYPE(20),IDK
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
C BK     1IDIMM
      DATA CHAR/' '/
C
C     CALCULATE PRESCRIBED T - O DISTANCES
C
      DO 35 NZ=1,NZA
      DO 15 I=1,4
      IF (ICON(NZ,I+1).EQ.0.AND.ICON(NZ,I+5).EQ.0) GO TO 20
      IF (ICON(NZ,I+1).EQ.0.AND.ICON(NZ,I+5).NE.0) GO TO 10
      NAB=MSYMAB(ICON(NZ,I+1))
      NATC=MSYMAB(ICON(NZ,1))
      M=M+1
      IF (M.GT.IDIMM) GO TO 100
      NN(M,1)=ICON(NZ,1)
      NN(M,2)=ICON(NZ,I+1)
      NN11(M,1)=NATC
      NN11(M,2)=NAB
      MM1=M-1
      IF (MM1.EQ.0) GO TO 2
      DO 1 J=1,MM1
      JS=J
      MS=M
      CALL EQUI
      IF (ITEST.EQ.1) GO TO 9
    1 CONTINUE
    2 IF (ICON(NZ,I+5).EQ.0) GO TO 13
      NAT=MSYMAB(ICON(NZ,I+5))
      NT1=ITYPE(NATC)
      NT2=ITYPE(NAB)
      NT3=ITYPE(NAT)
      IBTY=NT1+(NT2+2)*IDK+(NT3+2)*IDK*IDK
      DO 5 II=1,IA
      IF (BT(II).EQ.IBTY) GO TO 6
    5 CONTINUE
      WRITE (NTOUT,1015) KTYPE(NT1),KTYPE(NT2),KTYPE(NT3),NZ,ID(ICON(NZ,
     11)),ID(ICON(NZ,I+1)),ID(ICON(NZ,I+5))
      IERR=IERR+1
      GO TO 15
    6 IBT(NZ,I)=II
      DOB(M)=CAPID(II,1)
      DOBIN(M)=DOB(M)
      W(M)=WEIGHT(II,1)
      LJ(M)=M
      MW(M)=1
      ISYMB(M)=CHAR
      RN=RN+W(M)*W(M)*DOB(M)*DOB(M)
      MD(NZ,I)=M
      ICON(NZ,10)=ICON(NZ,10)+1
      GO TO 12
    9 M=M-1
      ICON(NZ,I+5)=-ICON(NZ,I+5)
      MD(NZ,I)=JS
      KK=I-1
      DO 11 K=1,KK
      IF (MD(NZ,K).EQ.JS) IBT(NZ,I)=IBT(NZ,K)
   11 CONTINUE
      GO TO 12
   10 WRITE (NTOUT,1020) NZ
      IERR=IERR+1
      GO TO 15
   13 WRITE (NTOUT,1035) NZ,ID(ICON(NZ,1)),ID(ICON(NZ,I+1))
      M=M-1
      IERR=IERR+1
   12 ICON(NZ,11)=ICON(NZ,11)+1
   15 CONTINUE
C
C     CALCULATE PRESCRIBED O - O DISTANCES
C
   20 K=4
      ITOT=ICON(NZ,11)
      IF (ICON(NZ,10).LT.2) GO TO 21
      IF (ITOT.EQ.4) GO TO 21
      WRITE (NTOUT,1025) NZ
      IERR=IERR+1
   21 IITO=ITOT+1
      DO 35 I=2,ITOT
      II=I+1
      DO 35 L=II,IITO
      K=K+1
      M=M+1
      IF (M.GT.IDIMM) GO TO 100
      NN(M,1)=ICON(NZ,I)
      NN(M,2)=ICON(NZ,L)
      NN11(M,1)=MSYMAB(ICON(NZ,I))
      NN11(M,2)=MSYMAB(ICON(NZ,L))
      MM1=M-1
      DO 25 J=1,MM1
      JS=J
      MS=M
      CALL EQUI
      IF (ITEST.EQ.1) GO TO 30
   25 CONTINUE
      DA=DOB(MD(NZ,I-1))
      DB=DOB(MD(NZ,L-1))
      DOB(M)=SQRT(DA*DA+DB*DB+2./3.*DA*DB)
      DOBIN(M)=DOB(M)
      W(M)=(WEIGHT(IBT(NZ,I-1),2)+WEIGHT(IBT(NZ,L-1),2))*.5
      LJ(M)=M
      MW(M)=1
      ISYMB(M)=CHAR
      RN=RN+W(M)*DOB(M)*DOB(M)
      MD(NZ,K)=M
      GO TO 35
   30 MD(NZ,K)=JS
      M=M-1
   35 CONTINUE
C
C     DETERMINE OUTER O - T DISTANCES
C
      DO 45 NZ=1,NZA
      ITOT=ICON(NZ,10)+1
      I=1
      DO 45 I1=2,ITOT
   37 I=I+1
      IF (ICON(NZ,I+4).LE.0) GO TO 37
      NAB=ICON(NZ,I)
      NAT=ICON(NZ,I+4)
      DO 40 MM=1,M
      IF (NN11(MM,1).EQ.MSYMAB(NAT).AND.NN11(MM,2).EQ.MSYMAB(NAB)) GO TO
     1 41
   40 CONTINUE
      WRITE (NTOUT,1030) ID(NAT),ID(NAB)
      IERR=IERR+1
      GO TO 45
   41 MD(NZ,I+9)=MM
   45 CONTINUE
C
C     CALCULATE PRESCRIBED T - T DISTANCES
C
      MM=M
      DO 60 NZ=1,NZA
      ITOT=ICON(NZ,10)+5
      I=5
      DO 60 I1=6,ITOT
   47 I=I+1
      IF (ICON(NZ,I).LE.0) GO TO 47
      NAT=ICON(NZ,I)
      M=M+1
      IF (M.GT.IDIMM) GO TO 100
      NN(M,1)=ICON(NZ,1)
      NN(M,2)=ICON(NZ,I)
      NN11(M,1)=MSYMAB(ICON(NZ,1))
      NN11(M,2)=MSYMAB(ICON(NZ,I))
      MM1=M-1
      DO 50 J=MM,MM1
      JS=J
      MS=M
      CALL EQUI
      IF (ITEST.EQ.1) GO TO 55
   50 CONTINUE
      DA=DOB(MD(NZ,I-5))
      DB=DOB(MD(NZ,I+5))
      CO=COS(CAPID(IBT(NZ,I-5),4)*PI180)
      DOB(M)=SQRT(DA*DA+DB*DB-2.*DA*DB*CO)
      DOBIN(M)=DOB(M)
      W(M)=WEIGHT(IBT(NZ,I-5),3)
      LJ(M)=M
      MW(M)=1
      ISYMB(M)=CHAR
      RN=RN+W(M)*DOB(M)*DOB(M)
      MD(NZ,I+9)=M
      GO TO 60
   55 MD(NZ,I+9)=JS
      M=M-1
   60 CONTINUE
      RETURN
  100 WRITE (NTOUT,1040) IDIMM
 1040 FORMAT (//'***','* DIMENSION CHECK : NUMBER LIMITED TO * ',I4)
      IERR=IERR+1
      M=M-1
      RETURN
 1015 FORMAT (//'***', 52H ERROR: NO DISTANCE FUNCTION SUPPLIED FOR BOND
     1TYPE ,A3,  3H - ,A3,  3H - ,A3/13X, 27HSPECIFIED ON TETCON CARD NO
     2,I3,  2H (,A6,  2H- ,A6,  4H... ,A6,  6H BOND))
 1020 FORMAT (//'*** ERROR ON TETCON CARD NO ',I3,
	1' : NUMBER OF BRIDGING ATOMS < NUMBER OF T - ATOMS')
 1025 FORMAT (//'***', 24H ERROR ON TETCON CARD NO,I3, 58H : ALL 4 O - AT
     1OMS ARE REQUIRED TO SPECIFY ALL INDEPENDENT, 16H O - O DISTANCES)
 1030 FORMAT (//'***', 50H ERROR: EQUIVALENT BOND IN ASYMMETRIC UNIT TO T
     1HE ,A6,  2H- ,A6, 39H BOND IS NOT SPECIFIED ON A TETCON CARD)
 1035 FORMAT (//'***', 24H ERROR ON TETCON CARD NO,I3,  3H : , 30HMISSING
     1 OUTER T-ATOM FOR BOND ,A6,  2H- ,A6)
      END
      SUBROUTINE RHO
      REAL*8 ID
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG,NEQU,
     1NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK3/ IBL1,ICAL,ICOR,IDIS,IDOB,IGIT,IHI,IHR,IMAT,INEW,IO
     1RT,IPAGE,IPCH,IZEIL,LETZTE,CORR
      COMMON /BLOCK4/ RN,RSIGMA,RWERT,ICVTST
      COMMON /BLOCK5/ D(400),DOB(400),DOBOLD(400),DELTAD(400),W(400),NN(
     1400,2),LJ(400),RAT(40),ISYMB(400),LDR(40)
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK7/ IDEC(100),KA(2,150),LCELL(6)
      COMMON /BLOCK8/ SI(3,150),SIGN(3,3,150),SYMOPS(10,200),B(3,3,200),
     1BI(3,200),IND,MTR(200)
      COMMON /BLOCK9/ DELV(300),IFELD(6),IV(300),IV1(300),KVAR(3,200)
      COMMON S(50000)
      COMMON /BLOC10/ T(300),G(3,3),GABL1(3,3,6),GABL2A(3,6,6),GABL2B(3,
     16,6),GABL2C(3,6,6),KRON(2,2),KRONEK(3,3),DGEL(300),MEIG
      COMMON /BLOC12/ SYMSIG(3,3,40),KATOR(400),NN11(400,2)
      COMMON /BLOC14/ IZ1(35),WR(35),BEDING(35),RESKO(35,20),KARI(35,20)
     1,LINH(35),LIND(35)
      COMMON /BLOC15/ NAPID,IRNG,KRNG
      COMMON /BLOC16/ ITYPE(200),BT(15),IA,NZA,ICON(40,11),MD(40,18),NCA
     1PID,CAPID(15,4),IAT,IAB,IBT(40,4),DOT,WEIGHT(15,3),DELTAX(3),DELTA
     20(3),DA,DB,ANG,TOT(40,4),OTO(40,6),PI180,IMOV(3),MSYMAB(200),DOBIN
     3(400),KTYPE(20),IDKTYP
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,MEQD,MCON,WFAC,
     1IDIMM,IDIMN,IDIMNZ,IDNEQU
C BK     1IDIMM
C BK      COMMON /BLOC18/ ROV,NC,NA,NS,IVF
      COMMON /BLOC18/ ROV,NC,NA,NS,IVF,IRNG1
      ROV=0.0
      DO 203 J=1,M
      IAB=NN(J,1)
      IAT=NN(J,2)
      CALL DIST
      D(J)=DOT
      JJ=LJ(J)
      DELTAD(J)=DOB(JJ)-D(J)
  203 ROV=ROV+W(J)*DELTAD(J)*DELTAD(J)
      IF (ICVTST.EQ.0.AND.NCYC.LT.NCYCLE) GO TO 201
      HN=SQRT(ROV)
      RWERT=HN/SQRT(RN)
      RSIGMA=HN/SQRT(FLOAT(M-NV))
  201 IF (NRE.EQ.NRH) RETURN
      ROL=0.0
      DO 204 L=1,NRE
      IF (LINH(L).EQ.1) GO TO 204
      ROL=ROL+WR(L)*BEDING(L)*BEDING(L)
  204 CONTINUE
      ROV=ROV+ROL
      RETURN
      END
      SUBROUTINE GENER
C***********************************************************************
C********* GENERATE SYMEQ CARDS  *********************************
C***********************************************************************
      REAL*8 ID
C BK B
      CHARACTER*8 CID
      DIMENSION CID(200)
      CHARACTER*4 CSYMOPS
      DIMENSION CSYMOPS(10,200)
C BK E
      REAL  IFAK,IGIN,IDEC,LCELL
      INTEGER BT
      CHARACTER*1 ITYP
      CHARACTER*10 KSYMB
      DIMENSION D(9)
      COMMON /BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
      COMMON /BLOCK2/ N,N1,N2,NV,NV1,NVD,NVG,NVO,NCYC,NCYCLE,NDIAG
     2 ,NEQU,NFUNC,NRE,M,MATCOR,NRAGA,NCC,NDO,NRH
      COMMON /BLOCK6/ ID(200),A(6),AOLD(6),X(3,200),XOLD(3,200),XINI(3,
     1200),FUDGE(3),IFAK(3,200),IGIN(6),ITIT(18),XMUL(200)
      COMMON /BLOCK7/ IDEC(100),KA(2,150),LCELL(6)
      COMMON /BLOCK8/ SI(3,150),SIGN(3,3,150),SYMOPS(10,200)
     1 ,B(3,3,200),BI(3,200),IND,MTR(200)
      COMMON /BLOC16/ ITYPE(200),BT(15), IA,NZA,ICON(40,11),
     $ MD(40,18),NCAPID,CAPID(15,4),IAT,IAB,IBT(40,4),DOT,WEIGHT(15,3),
     * DELTAXA(3),DELTAXB(3),DA,DB,ANG,TOT(40,4),OTO(40,6),PI180,IMOV(3)
     + ,MSYMAB(200),DOBIN(400),KTYPE(20),NATYP
      COMMON /BLOC17/ JS,MS,ITEST,MW(400),MEFF,WRF,IERR,
     $ MEQU,MCON,WFAC,IDIMM,IDIMN,IDIMNZ,IDNEQU
C BK      COMMON /BLOC22/ NUMB(10),IXYZ(3),LCENT(4),ISLASH,IBLANK,ICOMMA,
C BK     $ IPLUS,IMINUS,ISTAR
      COMMON /BLOC22/ IXYZ(3),LCENT(4),ISLASH,IBLANK,ICOMMA,
     $ IPLUS,IMINUS,ISTAR
      COMMON /BLOC30/ ISPACE,INPF,IPRINT,ICENT,NSYM,LTRAN,RL(3,4),
     $ KSYMB(3,48),RS(48,3,4),ITRAN(30),ISYMB(30)
      COMMON /BLOC31/ AXIS(3),KEY(3),DLIM(3),XI(3),XJS(3),
     1 T(3),SEP(3),ISGA(200),NAMX(6)
     2 ,RADI(20),ITET(200)
C BK B
      COMMON /BLOCKBK/ NUMB(10)
      EQUIVALENCE (CID,ID)
      EQUIVALENCE (CSYMOPS,SYMOPS)
C BK E
C     EQUIVALENCE (A(1),ACELL),(A(2),BCELL),(A(3),CCELL),
C    1 (A(4),ALF),(A(5),BET),(A(6),GAM)
C BK NUMB wird auch in SYMCRD initialisiert. Das mag f2c uerberhaupt nicht.
C BK Jetzt im Hauptprogramm:
C BK      DATA LBLANK,LSTAR,NUMB/1H ,1H*,1H1,1H2,1H3,1H4,1H5,1H6,1H7,
C BK     1 1H8,1H9,1H0/
      DATA LBLANK,LSTAR/1H ,1H*/
      DMINSQ = 1.E-10
      ACELL=A(1)
      BCELL=A(2)
      CCELL=A(3)
      ALF=A(4)
      BET=A(5)
      GAM=A(6)
      SA=SIN(ALF)
      SB=SIN(BET)
      SG=SIN(GAM)
      CA=COS(ALF)
      CB=COS(BET)
      CG=COS(GAM)
      CGSTAR= (CA*CB - CG)/(SA*SB)
      SGSTAR= SQRT(1.0 - CGSTAR**2)
      OR1=ACELL*SB*SGSTAR
      OR2= -ACELL*SB*CGSTAR
      OR3= BCELL*SA
      OR4=ACELL*CB
      OR5=BCELL*CA
      OR6= CCELL
      IF(IPRINT.EQ.0) GO TO 915
      PRINT 907,ACELL,BCELL,CCELL,ALF,BET,GAM
  907 FORMAT(//5X,' CELL PARAMETERS',6F10.3)
C BK      PRINT912
      PRINT 912
C BK 912   FORMAT(//' INPUT COORDS. FRACTIONAL '/)
  912 FORMAT(//' INPUT COORDS. FRACTIONAL '/)
      DO 910 I=1,N1
  910 PRINT 913, ID(I), (X(K,I),K=1,3)
  913 FORMAT(1H ,A6,3F10.4)
C
C**** SET UP LIMITS FOR TESTING IF ATOM IN RANGE FOR DISTANCE CALCULATIO
C
  915 VOL=ACELL*BCELL*CCELL*SQRT(1.0-CA*CA-CB*CB-CG*CG+2.0*CA*CB*CG)
      D100=VOL/(BCELL*CCELL*SA)
      D010=VOL/(ACELL*CCELL*SB)
      D001=VOL/(ACELL*BCELL*SG)
C**** DECIDE ORDER OF AXES IN DESCENDING AXIAL LENGTH.
      AXIS(1)=ACELL
      AXIS(2)=BCELL
      AXIS(3)=CCELL
      IF(ACELL.EQ.BCELL .OR. ACELL.EQ.CCELL)  AXIS(1)=ACELL+.0002
      IF(BCELL.EQ.CCELL) AXIS(2)=BCELL+.0001
      AXHI=AMAX1(AXIS(1),AXIS(2),AXIS(3))
      AXLO=AMIN1(AXIS(1),AXIS(2),AXIS(3))
      DO 70 K=1,3
      IF(AXIS(K).EQ.AXHI) GO TO 71
      IF(AXIS(K).EQ.AXLO)  GO TO 72
      KEY(K)=2
      GO TO 70
   71 KEY(K)=1
      GO TO 70
   72 KEY(K)=3
   70 CONTINUE
      KX=KEY(1)
      KY=KEY(2)
      KZ=KEY(3)
C
C     ASSIGN RADII TO DIFFERENT ATOM TYPES
C
      DO 210 NT = 1,NATYP
C BK      ENCODE(1,9200,ITYP) KTYPE(NT)
      WRITE(ITYP,9200) KTYPE(NT)
 9200 FORMAT(A1)
      IF(ITYP.EQ.'O')  GO TO 205
C     ASSIGN RADIUS OF 2.2 TO TETRAHEDRAL ATOMS (NON O ATOMS)
      RADI(NT) = 2.2
      GO TO 210
C     ASSIGN RADIUS OF 0.2 TO OXYGEN
  205 RADI(NT) = 0.2
  210 CONTINUE
C
C      TEST FOR TETRAHEDRAL ATOM
C
      DO 105 I=1,N1
      ISGA(I) = 0
      ITET(I) = 0
      IF(RADI(ITYPE(I)).LT.1.0) GO TO 105
      ITET(I) = 1
      NZA = NZA + 1
      IF(NZA.LE.IDIMNZ) GO TO 100
      WRITE(NTOUT,9230) IDIMNZ
 9230 FORMAT(//'***',' DIMENSION CHECK: NUMBER OF TETRAHEDRAL ATOMS '
     1'LIMITED TO'I4///'***',' PROGRAM STOPPED BY SUBROUTINE GENER '
     2 3H***)
      STOP
  100 ICON(NZA,1) = I
      DO 104 II=2,11
  104 ICON(NZA,II) = 0
  105 CONTINUE
      NAT1 = N1 + 1
C
C     INPUT SYMMETRY CARDS FORM FILE NTIN
C
      CALL SYMCRD
C
      IF(IPRINT.NE.0)  PRINT 9300
C
C**** START LOOP ON TETRAHEDRAL ATOM I
      DO 306 NZ = 1,NZA
      I = ICON(NZ,1)
      XI(KX)=X(1,I)
      XI(KY)=X(2,I)
      XI(KZ)=X(3,I)
 9300 FORMAT(//  5X,'ATOM I  ATOM J   DIJ    ',12X,'SYM TX TY TZ' ,12X,
     1  '   XI      YI      ZI       XJST    YJST    ZJST' /)
      LINE=LINE+2
C**** START LOOP ON ATOMS J
      DO 305 J=1,N1
      ISEQ = 0
      DMAX=RADI(ITYPE(I))+RADI(ITYPE(J))+0.25
      DMAXSQ=DMAX*DMAX
C**** SET UP THE D LIMITS FOR EACH AXIS
C
      DLIM(KX)=DMAX/D100
      DLIM(KY)=DMAX/D010
      DLIM(KZ)=DMAX/D001
C**** START LOOP ON SYMMETRY OPERATORS
      DO 304 NS= 1,NSYM
      NSS = NS
C**** GET XJS,YJS,ZJS
      XX=X(1,J)
      YY=X(2,J)
      ZZ=X(3,J)
      XF = RS(NS,1,1)*XX + RS(NS,1,2)*YY + RS(NS,1,3)*ZZ + RS(NS,1,4)
      YF = RS(NS,2,1)*XX + RS(NS,2,2)*YY + RS(NS,2,3)*ZZ + RS(NS,2,4)
      ZF = RS(NS,3,1)*XX + RS(NS,3,2)*YY + RS(NS,3,3)*ZZ + RS(NS,3,4)
      DO 304 IC = 1,ICENT
      IF(IC.EQ.1) GO TO 30
      XF = -XF
      YF = -YF
      ZF = -ZF
      NSS = - NSS
C***********************************************************************
C**** BEGIN TRANSLATIONAL SEARCH ROUTINE.
C***********************************************************************
   30 DO 304 LT = 1,LTRAN
      XJS(KX) = XF
      XJS(KY) = YF
      XJS(KZ) = ZF
      T(KX) = RL(1,LT)
      T(KY) = RL(2,LT)
      T(KZ) = RL(3,LT)
      DO 35 K = 1,3
   35 XJS(K) = XJS(K) + T(K)
      K=1
      IDIST = 0
C**** TRANSLATE XJS ON AXIS K TO A POINT BELOW LOWER LIMIT.
   27 IF(XI(K)-XJS(K) .GT.DLIM(K))  GO TO 20
      XJS(K)=XJS(K)-1.0
      T(K)=T(K)-1.0
      GO TO 27
C**** ADVANCE +1.0 ON AXIS K.
   20 XJS(K)=XJS(K)+1.0
      T(K)=T(K)+1.0
      SEP(K)=XJS(K) - XI(K)
      IF(SEP(K).LT. -DLIM(K)) GO TO 20
C**** TEST IF XJS-XI LESS THAN OR EQUAL TO DLIM(K)
      IF(SEP(K).GT.DLIM(K))  GO TO 23
      GO TO 24
C**** TEST GIVES ATOM OUT OF RANGE. REDUCE AXIS K BY 1.
   23 IF(IDIST.EQ.0)  GO TO 50
      K=K-1
      IF(K) 50,50,20
C**** K=0 SIGNALS THE END OF THE SEARCH AND WE LEAVE THROUGH 50.
C**** IF XJS WAS WITHIN RANGE ON AXIS K THEN WE ADVANCE K BY 1
   24 K=K+1
      IF(K.LE.3)  GO TO 27
C**** ATOM IN RANGE FOR ALL 3 AXES. NOW CALCULATE EXACT DISTANCE.
C**** GET DIJ.
   29 DXO=OR1*SEP(KX)
      DYO=OR2*SEP(KX) + OR3*SEP(KY)
      DZO=OR4*SEP(KX) + OR5*SEP(KY) + OR6*SEP(KZ)
      DIJSQ=DXO*DXO + DYO*DYO + DZO*DZO
      IDIST = 1
      IF(DIJSQ.LE.DMINSQ)  GO TO 32
      IF(DIJSQ.LE.DMAXSQ)  GO TO 360
   32 K=K-1
      GO TO 20
C**** WE HAVE A DISTANCE LESS THA OR EQUAL TO DMAX. PRINT AND STORE.
  360 DIJ=SQRT(DIJSQ)
      XJST=XJS(KX)
      YJST=XJS(KY)
      ZJST=XJS(KZ)
      IF(IPRINT.NE.0) PRINT 9360,ID(I),ID(J),DIJ,NSS,T(KX),T(KY),
     1T(KZ),XI(KX),XI(KY),XI(KZ),XJST,YJST,ZJST
      IF(DIJ.GT.0.05)  GO TO 361
      IF(IPRINT.EQ.0) GO TO 32
      IF(I.NE.J) PRINT 9366, ID(I),ID(J)
      IF(I.EQ.J)  PRINT 9367, ID(I)
 9360 FORMAT(5X,A6,2X,A6,F10.5,10X,I5,3F4.1,7X,2(2X,3F8.5))
 9366 FORMAT(' NOTE****',5X,A6,2X,A6,' ARE EQUIVALENT BY SPACE GRO
     1UP SYMMETRY')
 9367 FORMAT(' NOTE****',5X,A6,' LIES ON A SPECIAL POSITION')
      GO TO 32
C**** TEST FOR SYMMETRY GENERATED ATOM J'
  361 IF(ABS(X(1,J) -XJST) .GT. 0.0001)  GO TO 363
      IF(ABS(X(2,J)-YJST) .GT. 0.0001)  GO TO 363
      IF(ABS(X(3,J)-ZJST) .GT. 0.0001)  GO TO 363
      IF(ISEQ.EQ.1)  GO TO 32
      ISEQ = 1
      NJ = J
      GO TO 362
  363 IF(NAT1.GT.N) GO TO 280
      DO 250 L=NAT1,N
      IF(ABS(X(1,L)-XJST).GT.1.E-4) GO TO 250
      IF(ABS(X(2,L)-YJST).GT.1.E-4) GO TO 250
      IF(ABS(X(3,L)-ZJST).GT.1.E-4) GO TO 250
      NJ = L
      GO TO 290
  250 CONTINUE
C
C     STORE NEW SYMMETRY GENERATED ATOM
C
  280 N2 = N2 + 1
      N = N1 + N2
      IF(N.LE.IDIMN)  GO TO 282
      WRITE(NTOUT,9205) IDIMN
 9205 FORMAT(//'***',' DIMENSION CHECK: TOTAL NUMBER OF ATOMS LIMITED'
     1 ' TO'I4)
      IERR = IERR + 1
      N2 = N2 - 1
      N = N - 1
      GO TO 304
  282 NEQU = NEQU + 1
      IF (NEQU.LE.IDNEQU)  GO TO 284
      WRITE(NTOUT,9210) IDNEQU
 9210 FORMAT(//'***',' DIMENSION CHECK: NUMBER OF ATOMS GENERATED'
     1 ' EXEEDS LIMIT OF'I4)
      IERR = IERR + 1
      NEQU = NEQU - 1
      GO TO 304
  284 KA(1,NEQU) = J
      ISGA(J) = ISGA(J) + 1
      ISGA(N) = ISGA(J)
      KA(2,NEQU) = N
      MSYMAB(N) = KA(1,NEQU)
      MTR(N) = N
C BK B Unter Ausnutzung der Tatsache, dass (SYMOPS(I,N),I=6,10) mit
C      Blanks gefuellt wird:
C       ENCODE(50,9215,SYMOPS(1,N)) NSS,T(KX),T(KY),T(KZ)
      WRITE(CSYMOPS(1,N),9216) NSS
      DO 9217 IBK=2,5
       CSYMOPS(IBK,N) = ' '
 9217 CONTINUE
 9215 FORMAT(I4,16X,3(F5.2,5X))
 9216 FORMAT(I4)
C BK E
      DO 286 IS = 6,10
C BK  286 SYMOPS(IS,N) = 1H
  286 CSYMOPS(IS,N) = ' '
      NJ = N
      X(1,N) = XJST
      X(2,N) = YJST
      X(3,N) = ZJST
      SK = 1.0
      IF(IC.EQ.2 )  SK = -1.
      DO 287 M1 = 1,3
      DO 287  M2 = 1,3
  287 SIGN(M1,M2,NEQU) = SK * RS(NS,M1,M2)
      SI(1,NEQU) = SK * RS(NS,1,4) + T(KX)
      SI(2,NEQU) = SK * RS(NS,2,4) + T(KY)
      SI(3,NEQU) = SK * RS(NS,3,4) + T(KZ)
      DO 289 M1 = 1,3
      DO 289 M3 = 1,3
      BIK = 0.0
      DO 288 M2=1,3
  288 BIK = BIK + SIGN(M1,M2,NEQU) * B(M2,M3,J)
  289 B(M1,M3,N) = BIK
  290 IF(IPRINT.NE.0) PRINT 9361, ID(J)
 9361 FORMAT(' NOTE**** CONNECTION TO A SYMMETRY RELATED ATOM (',A6
     1'*')
C
C SET UP CONNECTION TABLE ICON AND NCON
C
C
  362 IF(ITET(J).EQ.0)  GO TO 120
      JT = ICON(I,11) + 5
      IF(JT.EQ.5)  GO TO 110
      DO 106 JJ = 6,JT
      IF(ICON(NZ,JJ).EQ.NJ)  GO TO 32
  106 CONTINUE
  110 JT = JT + 1
      IF(JT.GT.9)  GO TO 115
      ICON(NZ,11) = JT - 5
      ICON(NZ,JT) = NJ
      D(JT) = DIJSQ
      GO TO 32
  115 JTA = 6
      JTE = 9
      ICON(NZ,11) = 4
      GO TO 140
  120 JO = ICON(NZ,10) + 1
      IF(JO.EQ.1)  GO TO 130
      DO 125 JJ = 2,JO
      IF(ICON(NZ,JJ).EQ.NJ)  GO TO 32
  125 CONTINUE
  130 JO = JO + 1
      IF(JO.GT.5)  GO TO 135
      ICON(NZ,10) = JO -1
      ICON(NZ,JO) = NJ
      D(JO) = DIJSQ
      GO TO 32
  135 JTA = 2
      JTE = 5
      ICON(NZ,10) = 4
  140 DO 150 JT = JTA,JTE
      IF(D(JT).LT.DIJSQ)  GO TO 150
      D(1) = D(JT)
      D(JT) = DIJSQ
      DIJSQ = D(1)
      NJSAVE = ICON(NZ,JT)
      ICON(NZ,JT) = NJ
      NJ = NJSAVE
  150 CONTINUE
      GO TO 32
  160 WRITE(NTOUT,9260) ID(I),ID(J)
      IERR = IERR + 1
      GO TO 32
 9260 FORMAT(//'***',' ERROR IN SETTING UP OF CONNECTION TABLE'
     1 ' FOR ATOMS 'A6' AND 'A6)
  170 GO TO 32
C**** END OF TRANSLATIONAL SEARCH.......................................
   50 CONTINUE
C**** END LOOP ON SYMMETRY OPRS.
  304 CONTINUE
C**** END OF LOOP ON ATOMS J
  305 CONTINUE
  306 CONTINUE
C**** END OF LOOP ON I
C**** PRINT OUT LIST OF SYMMETRY GENERATED CONNECTIONS IF PRESENT.
      IF(NEQU.EQ.0)  GO TO 395
      IF(IPRINT.NE.0) PRINT 9393
 9393 FORMAT(//' LIST OF GENERATED SYMMETRY EQUIVALENT ATOMS'/)
      NAT = N1
      DO 394 L=1,NEQU
      J = KA(1,L)
      NAT = NAT + 1
C BK      DECODE(6,9395,ID(J))(NAMX(K),K=1,6)
      READ(CID(J),9395)(NAMX(K),K=1,6)
 9395 FORMAT(6A1)
      DO 391 K=2,6
      IF(NAMX(K).EQ.IBLANK) GO TO 392
  391 CONTINUE
      K=6
  392 NSTAR=ISGA(NAT)
      DO 425 NR=1,NSTAR
      IF(K.GT.6) K=6
      IF(NAMX(K).NE.IBLANK) GO TO 399
      NAMX(K)=ISTAR
      GO TO 420
  399 IF(NAMX(K).EQ.ISTAR) GO TO 410
      DO 400 NU=2,10
      IF(NAMX(K).EQ.NUMB(NU)) GO TO 405
  400 CONTINUE
      NU = 9
  405 NU = NU + 1
      GO TO 415
  410 NU = 2
  415 NAMX(K) = NUMB(NU)
  420 K=K+1
  425 CONTINUE
      IF(IPRINT.NE.0) PRINT 9394,NAMX,(SYMOPS(K,NAT),K=1,5),(X(K,NAT),
     1 K=1,3)
C BK      ENCODE(6,9395,ID(NAT))(NAMX(K),K=1,6)
      WRITE(CID(NAT),9395)(NAMX(K),K=1,6)
  394 CONTINUE
 9394 FORMAT(10X,6A1,5A10,5X,3F8.4)
  395 CONTINUE
      DMAXSQ = 7.0
      IF(IPRINT.NE.0) PRINT 9500
 9500 FORMAT(//' CONNECTIVITY ARRAYS (AS ON TETCON CARDS)'/)
      DO 550 NZ = 1,NZA
      IO = ICON(NZ,10)
      IF(IO.NE.4)  GO TO 600
      NT = ICON(NZ,11)
      IF(NT.NE.4)  GO TO 600
      DO 540 I = 1,IO
      NO = ICON(NZ,I+1)
  530 NTP1 = I + 5
      DMINSQ = DMAXSQ
      DO 535 K = NTP1,9
      K1 = ICON(NZ,K)
      SEP(1) = X(1,K1) - X(1,NO)
      SEP(2) = X(2,K1) - X(2,NO)
      SEP(3) = X(3,K1) - X(3,NO)
      DXO = OR1*SEP(1)
      DYO = OR2*SEP(1) + OR3*SEP(2)
      DZO = OR4*SEP(1) + OR5*SEP(2) + OR6*SEP(3)
      DIJSQ = DXO*DXO + DYO*DYO + DZO*DZO
      IF(DIJSQ.GT.DMINSQ)  GO TO 535
      DMINSQ = DIJSQ
      NSAVE = ICON(NZ,NTP1)
      ICON(NZ,NTP1) = K1
      ICON(NZ,K) = NSAVE
  535 CONTINUE
  540 CONTINUE
      IF(IPRINT.NE.0) WRITE(NTOUT,9550) (ID(ICON(NZ,K)),K=1,9),
     $ ICON(NZ,10),ICON(NZ,11)
      ICON(NZ,10) = 0
      ICON(NZ,11) = 0
 9550 FORMAT(1X,9A7,2I5)
      GO TO 550
  600 WRITE(NTOUT,9600)  ID(ICON(NZ,1)),(ID(ICON(NZ,IK)),IK=1,9),
     $ ICON(NZ,10),ICON(NZ,11)
      IERR = IERR + 1
  550 CONTINUE
      IF(IPRINT.NE.0)  CALL KOPF
      RETURN
 9600 FORMAT(//'***',' ERROR IN SETTING UP OF CONNECTION TABLE FOR '
     1 'ATOM 'A6,/' ICON - ARRAY = ',9A7,2I5)
      END
      SUBROUTINE SYMCRD
      CHARACTER*10 KSYMB
      CHARACTER*30 KXSYM
      COMMON / BLOCK1/ NTIN,NTOUT,NFILEA,NFILEB,IFLAG
C BK      COMMON /BLOC22/ NUMB(10),IXYZ(3),LCENT(4),ISLASH,IBLANK,ICOMMA,
C BK     $ IPLUS,IMINUS,ISTAR
      COMMON /BLOC22/ IXYZ(3),LCENT(4),ISLASH,IBLANK,ICOMMA,
     $ IPLUS,IMINUS,ISTAR
      COMMON /BLOC30/ ISPACE,INPF,IPRINT,ICENT,NSYM,LTRAN,RL(3,4),
     $ KSYMB(3,48),RS(48,3,4),ITRAN(30),ISYMB(30)
      COMMON /BLOCKBK/ NUMB(10)
C BK wird schon im Hauptprogramm initialisiert, macht aber nichts:
      DATA IXYZ/1HX,1HY,1HZ/
      DATA ISTAR/1H*/
C BK Wird auch in GENER initialisiert. Das mag f2c uerberhaupt nicht.
C BK Jetzt im Hauptprogramm:
C BK      DATA NUMB/1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H0/
      DATA ISLASH,IBLANK,ICOMMA,IPLUS,IMINUS/1H/,1H ,1H,,1H+,1H-/
      DATA LCENT/'ACEN','TRIC','CENT','RIC '/
      REWIND NTIN
    5 READ(NTIN,100) ISG,NSYM,ICENT,LTRAN,(ITRAN(I),I=1,30)
  100 FORMAT(I5,3I3,30I1)
      IF(ISG.EQ.ISPACE) GO TO 20
      DO  10 I = 1,NSYM
   10 READ(NTIN,110) IDUM
  110 FORMAT(I4)
      GO TO 5
   20 ICENT = ICENT + 1
      LTRAN = LTRAN + 1
      I = 1
      DO 23 K=1,3
   23 RL(K,1) = 0.0
      IF(LTRAN.EQ.1)   GO TO 26
      DO 25 L = 2,LTRAN
      DO 25 K = 1,3
      IF(ITRAN(I+1).NE.0) GO TO 24
      RL(K,L) = 0.0
      GO TO 25
   24 RL(K,L) = FLOAT(ITRAN(I))/FLOAT(ITRAN(I+1))
   25 I = I + 2
   26 DO 35 N = 1,NSYM
      READ(NTIN,120) (ITRAN(M),M=1,16)
  120 FORMAT(3(I4,2I2,2I3,1X),I3)
      IF(ITRAN(16).NE.N)  GO TO 900
      M = 1
      DO 30 K=1,3
      IF(ITRAN(M+1).NE.0) GO TO 27
      RS(N,K,4) = 0.0
      GO TO 28
   27 RS(N,K,4) = FLOAT(ITRAN(M))/FLOAT(ITRAN(M+1))
   28 M = M + 2
      DO 30 L = 1,3
      RS(N,K,L) = ITRAN(M)
      M = M + 1
   30 CONTINUE
      DO 310 KK=1,30
  310 ISYMB(KK) = IBLANK
      K = 0
      I = 1
      DO 360 IS = 1,3
      IF(ITRAN(I).EQ.0) GO TO 320
      K = K+1
      ISYMB(K) = NUMB(ITRAN(I))
      K = K+1
      ISYMB(K) = ISLASH
      K = K + 1
      ISYMB(K) = NUMB(ITRAN(I+1))
  320 I = I + 2
      DO 350 KX =1,3
      IF(ITRAN(I).EQ.0)  GO TO 350
      IF(ITRAN(I).GE.0)  GO TO 330
  325 K = K + 1
      ISYMB(K) = IMINUS
      ITRAN(I) = -ITRAN(I)
      K = K + 1
      GO TO 335
  330 K = K + 1
      IF(K.EQ.1)  GO TO 335
      IF(ISYMB(K-1).EQ.IBLANK)  GO TO 335
      ISYMB(K) = IPLUS
      K = K + 1
  335 ISYMB(K) = IXYZ(KX)
  350 I = I+1
      IF(IS.EQ.3)  GO TO 360
      K = IS*10
      ISYMB(K) = ICOMMA
      K = K + 1
      ISYMB(K) = IBLANK
  360 CONTINUE
C BK B    ENCODE(30,370,KSYMB(1,N)) (ISYMB(KK),KK=1,K)

      IF (K.GT.20) THEN
        WRITE(KSYMB(3,N),370) (ISYMB(KK),KK=21,K)
      ENDIF

      IF (K.GT.10) THEN
        WRITE(KSYMB(2,N),370) (ISYMB(KK),KK=11,MIN(K,20))
      ENDIF

      WRITE(KSYMB(1,N),370) (ISYMB(KK),KK=1,MIN(K,10))

C   370 FORMAT(30A1)
  370 FORMAT(10A1)
C BK E
   35 CONTINUE
      IF(IPRINT.EQ.0)  RETURN
      WRITE(NTOUT,190) ISPACE,NTIN
  190 FORMAT(/'READ SYMMETRY INFORMATION FOR SPACE GROUP NR'I5,
     1 ' FROM FILE'I3/)
      IF(ICENT.EQ.2) ICENT=3
      WRITE(NTOUT,200)  ISPACE,(LCENT(I),I=ICENT,ICENT+1)
  200 FORMAT(//'SPACE GROUP NR 'I5' IS '2A4)
      IF(LTRAN.EQ.1)  GO TO 130
      WRITE(NTOUT,205) ((RL(K,L),K=1,3),L=1,LTRAN)
  205 FORMAT(/'LATTICE TRANSLATIONS'///4(3F4.1,10X))
  130 WRITE(NTOUT,212) (N,(KSYMB(I,N),I=1,3),N=1,NSYM)
  212 FORMAT(///'CODED SYMMETRY TRANSFORMATIONS'//(1X,I5,3X,3A10))
      RETURN
  900 WRITE(NTOUT,910) NTIN,ISPACE
  910 FORMAT(//'***',' ERROR IN READING THE SYMMETRY INFORMATION'
     1 ' FROM FILE'I5' FOR SPACE GROUP NR'I6)
      STOP
      END

      block data ranini

      real U(97), C, CD, CM
      integer I97, J97
      logical TEST
      common /raset1/ U, C, CD, CM, I97, J97, TEST

      data TEST /.FALSE./

      end


C This random number generator originally appeared in "Toward a Universal
C Random Number Generator" by George Marsaglia and Arif Zaman.
C Florida State University Report: FSU-SCRI-87-50 (1987)
C
C It was later modified by F. James and published in "A Review of Pseudo-
C random Number Generators"
C
C THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
C       (However, a newly discovered technique can yield
C         a period of 10^600. But that is still in the development stage.)
C
C It passes ALL of the tests for random number generators and has a period
C   of 2^144, is completely portable (gives bit identical results on all
C   machines with at least 24-bit mantissas in the floating point
C   representation).
C
C The algorithm is a combination of a Fibonacci sequence (with lags of 97
C   and 33, and operation "subtraction plus one, modulo one") and an
C   "arithmetic sequence" (using subtraction).
C
C On a Vax 11/780, this random number generator can produce a number in
C    13 microseconds.
C========================================================================

      subroutine RMARIN(IJ,KL)
C This is the initialization routine for the random number generator RANMAR()
C NOTE: The seed variables can have values between:    0 <= IJ <= 31328
C                                                      0 <= KL <= 30081
C The random number sequences created by these two seeds are of sufficient
C length to complete an entire calculation with. For example, if sveral
C different groups are working on different parts of the same calculation,
C each group could be assigned its own IJ seed. This would leave each group
C with 30000 choices for the second seed. That is to say, this random
C number generator can create 900 million different subsequences -- with
C each subsequence having a length of approximately 10^30.
C
C Use IJ = 1802 & KL = 9373 to test the random number generator. The
C subroutine RANMAR should be used to generate 20000 random numbers.
C Then display the next six random numbers generated multiplied by 4096*4096
C If the random number generator is working properly, the random numbers
C should be:
C           6533892.0  14220222.0  7275067.0
C           6172232.0  8354498.0   10633180.0


      real U(97), C, CD, CM
      integer I97, J97
      logical TEST
      common /raset1/ U, C, CD, CM, I97, J97, TEST

      if( IJ .lt. 0  .or.  IJ .gt. 31328  .or.
     *    KL .lt. 0  .or.  KL .gt. 30081 ) then
          print '(A)', ' The first random number seed must have a value
     *between 0 and 31328'
          print '(A)',' The second seed must have a value between 0 and     
     *30081'
            stop
      endif

      i = mod(IJ/177, 177) + 2
      j = mod(IJ    , 177) + 2
      k = mod(KL/169, 178) + 1
      l = mod(KL,     169)

      do 2 ii = 1, 97
         s = 0.0
         t = 0.5
         do 3 jj = 1, 24
            m = mod(mod(i*j, 179)*k, 179)
            i = j
            j = k
            k = m
            l = mod(53*l+1, 169)
            if (mod(l*m, 64) .ge. 32) then
               s = s + t
            endif
            t = 0.5 * t
C BK 3        continue
 3       continue
         U(ii) = s
C BK 2     continue
 2    continue

      C = 362436.0 / 16777216.0
      CD = 7654321.0 / 16777216.0
      CM = 16777213.0 /16777216.0

      I97 = 97
      J97 = 33

      TEST = .TRUE.
      return
      end

      subroutine ranmar(RVEC, LEN)
C This is the random number generator proposed by George Marsaglia in
C Florida State University Report: FSU-SCRI-87-50
C It was slightly modified by F. James to produce an array of pseudorandom
C numbers.
      REAL RVEC(*)
      real U(97), C, CD, CM
      integer I97, J97
      logical TEST
      common /raset1/ U, C, CD, CM, I97, J97, TEST

      integer ivec

      if( .NOT. TEST ) then
         print '(A)',' Call the init routine (RMARIN) before calling RAN    
     *MAR'
         stop
      endif

      do 100 ivec = 1, LEN
         uni = U(I97) - U(J97)
         if( uni .lt. 0.0 ) uni = uni + 1.0
         U(I97) = uni
         I97 = I97 - 1
         if(I97 .eq. 0) I97 = 97
         J97 = J97 - 1
         if(J97 .eq. 0) J97 = 97
         C = C - CD
         if( C .lt. 0.0 ) C = C + CM
         uni = uni - C
         if( uni .lt. 0.0 ) uni = uni + 1.0
         RVEC(ivec) = uni
C BK 100   continue
 100  continue
      return
      end

