C 19 SEP 05 - MWS - ADD TRUE NUCLEAR CHARGE ARRAY TO INFOA COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 30 APR 05 - DGF - SYNCHRONISE INFO COMMON BLOCK
C  7 NOV 00 - MWS - MBLDR: TRAP IF MXATM IS EXCEEDED
C 21 DEC 99 - MWS - MBLDR: READAT DUMMY ARG, REWRITE PRAXIS ROTATION
C 29 AUG 99 - JPG - ZMATAB: FIX BUG WITH NEGATIVE ANGLES
C 12 NOV 98 - GDF - MAKE SYMTRY COMMON CONSISTENT
C  5 MAR 96 - MWS - CHANGE SHELL SYMMETRY PACKING COMMONS
C  1 FEB 95 - MWS - ZMATAB: INITIALIZE EOF BEFORE READING
C 11 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 27 OCT 94 - MWS - PRAXIS: PAD OUT MACHSW COMMON
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C  1 JUN 94 - MWS - MBLDR: LET COORD=CART PROCEED IF NZVAR.GT.0
C 13 APR 94 - MWS - RDMPC,RDZMT: COUNT Z-MATRIX VALUES INPUT
C  3 NOV 93 - MH  - MBLDR: CHARMM READS ATOMS BY A SPECIAL CALL
C  1 NOV 93 - MWS - PRAXIS: USE PROPER ROTATION MATRIX
C 29 SEP 93 - MWS - RDMPC,ZMATAB: CHECK INPUT ATOM CONNECTIVITY
C 16 JUN 93 - MWS - RDMPC: DEFAULT CONNECTION OF 3RD ATOM IS TO 2ND
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C 16 FEB 92 - JHJ - RDMPC: CONNECTIVITY FOR ATOM 2 READ IN.
C 30 JAN 92 - TLW - RDZMT: FIX INTERNAL READS
C 26 JAN 92 - MWS - EQVATM: DELETE AUTOMATIC Z-MATRIX TEST
C 11 JAN 92 - TLW - RDZMT: MAKE READS PARALLEL
C 10 JAN 92 - TLW - MBLDR,RDMPC,ZMATAB: CALL ABRT BEFORE STOP
C  7 JAN 92 - TLW - MAKE WRITES PARALLEL;ADD COMMON PAR
C  5 DEC 91 - MWS - RDZMT: PRINT ERROR MESSAGE BEFORE ABRT
C 28 FEB 91 - MWS - PRAXIS: HARDWIRE JACOBI DIAGAONLIZATION
C  9 NOV 90 - JHJ - NEW MODULE ADDED FOR Z-MATRIX INPUT IN $DATA
C
C*MODULE INPUTC   *DECK EQVATM
      SUBROUTINE EQVATM(AZNUC,CORD,CORD2,IANS,NAT,I1,I2,I3)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, TM4=1.0D-04, TM12=1.0D-12)
C
      DIMENSION CORD(MXATM,3),CORD2(MXATM,3),AZNUC(MXATM)
C
C     ----- TEST FOR IDENTICAL MOLECULES -----
C     THE TWO SETS OF COORDINATES ARE CORD,CORD2.
C     THE PERMUTATION BEING CONSIDERED IS I1,I2,I3.
C     RETURN IANS=0 IF BOTH SETS MATCH EXACTLY.
C
      DO 100 I = 1,NAT
         DO 50 J = 1,NAT
            DZN = (AZNUC(I) - AZNUC(J))**2
            IF(DZN.GT.TM4) GO TO 50
            TEST=(CORD(I,I1)-CORD2(J,1))**2
     *          +(CORD(I,I2)-CORD2(J,2))**2
     *          +(CORD(I,I3)-CORD2(J,3))**2
            IF(TEST.LE.TM12) GO TO 100
   50    CONTINUE
         IANS=1
         RETURN
  100 CONTINUE
      IANS = 0
      RETURN
      END
C*MODULE INPUTC   *DECK MBLDR
      SUBROUTINE MBLDR(IUNTRD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION IXX,IYY,IZZ
C
      LOGICAL GAUSS, MOPAC, INTRNL
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXSH=5000)
      PARAMETER (MXCHRM=1)
C
      CHARACTER*10 ENDWRD, AATOM, UATOM
C
      COMMON /CHMGMS/ XCHM(MXCHRM),YCHM(MXCHRM),ZCHM(MXCHRM),
     *                DXELMM(MXCHRM),DYELMM(MXCHRM),DZELMM(MXCHRM),
     *                QCHM(MXCHRM),NCHMAT,KCHRMM
      COMMON /COORDN/ UZNUC(MXATM),CUNIQ(MXATM,3),NATREL,UATOM(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TRANSF/ XOLD,YOLD,ZOLD,XNEW,YNEW,ZNEW,XP,YP,ZP
      COMMON /ZMTALT/ NZMAT2,NZVAR2,NVAR2,NZMTRD,ICOORD
      COMMON /ZMINF / CORD(MXATM,3),AZNUC(MXATM),CORD2(MXATM,3),
     *                AATOM(MXATM)
C
      PARAMETER (TM8=1.0D-08)
C
      DATA ENDWRD /'$END      '/
C
C       ICOORD=1,2,3 MEANS COORD=CART,ZMT,ZMTMPC
C       IUNTRD=-1,1 MEANS UNITS=BOHR,ANGS
C       THESE ARE THE ONLY POSSIBLE VALUES ON ENTRY.
C
      IANS = 0
      NAT = 0
C
      GAUSS = ICOORD.EQ.2
      MOPAC = ICOORD.EQ.3
      INTRNL = GAUSS  .OR.  MOPAC
C
C         READ IN COORDINATES IN Z-MATRIX FORM
C
      IF(INTRNL) THEN
         NZVAR2 = 0
         CALL ZMATAB(MOPAC,NAT)
         GO TO 200
      END IF
C
C        READ ALL CARTESIAN COORDINATES
C        DMY ARGUMENT IS OK SINCE THIS CALL NEVER TRIES FOR HINT COORDS
C
      IF(KCHRMM.EQ.0) THEN
  110    CONTINUE
            NAT = NAT + 1
            IF(NAT.GT.MXATM) GO TO 810
            DMY = 0.0D+00
            CALL READAT(AATOM(NAT),AZNUC(NAT),CORD(NAT,1),CORD(NAT,2),
     *                  CORD(NAT,3),DMY,NAT,ICOORD,IUNTRD)
            IF(AATOM(NAT).EQ.ENDWRD) GO TO 120
         GO TO 110
C
  120    CONTINUE
         NAT=NAT - 1
      ELSE
         CALL CHMDAT(AATOM,AZNUC,CORD,NAT)
         RETURN
      END IF
C
C     ----- DISCARD ANY DUMMY ATOMS -----
C     IN GAMESS DUMMY ATOMS MUST APPEAR LAST, AND SO CAN ONLY BE
C     USED WITH COORD=UNIQUE AND HINT RUNS, WHERE THEY MUST BE
C     INPUT *AT THE END* OF ANY REAL ATOMS.   HERE IT IS SIMPLER
C     JUST TO OMIT THE "STUPID" THINGS.
C
  200 CONTINUE
      NATREL = 0
      DO 210 I=1,NAT
         NUCZ = INT(AZNUC(I))
         IF(NUCZ.GT.0) THEN
            NATREL = NATREL+1
            UATOM(NATREL) = AATOM(I)
            UZNUC(NATREL) = AZNUC(I)
            CUNIQ(NATREL,1) = CORD(I,1)
            CUNIQ(NATREL,2) = CORD(I,2)
            CUNIQ(NATREL,3) = CORD(I,3)
         END IF
  210 CONTINUE
C
      IF(NATREL.LT.NAT) NATREL = -NATREL
      NAT = IABS(NATREL)
C
      DO 220 I=1,NAT
         AATOM(I) = UATOM(I)
         AZNUC(I) = UZNUC(I)
         CORD(I,1) = CUNIQ(I,1)
         CORD(I,2) = CUNIQ(I,2)
         CORD(I,3) = CUNIQ(I,3)
  220 CONTINUE
C
C     ----- COMPUTE PRINCIPAL AXES -----
C     "TALKING OF AXES,"
C     SAID THE DUTCHESS,
C     "OFF WITH HER HEAD!"
C        -LEWIS CARROLL
C
      CALL PRAXIS(CORD,AZNUC,NAT,MXATM,IXX,IYY,IZZ)
      IF (MASWRK) WRITE(IW,9020) IXX,IYY,IZZ
C
C SWITCH AXES IF NECESSARY; LIKE SO..
C 1.ALL SIX COMBINATIONS OF X,Y,Z ARE,IN TURN, SUBJECTED TO ALL
C   NT SYMMETRY OPERATIONS.
C 2.AND COMPARED TO THE INITIAL STRUCTURE.
C 3.IF INITIAL STRUCTURE.EQ.FINAL STRUCTURE AXES ARE CORRECT.
C WE'RE OPTIMISTIC SO X,Y,Z IS TRIED FIRST.
C
      DO 430 IAX = 1,3
         IICNT = 0
         I1 = IAX
         I2 = 1 + MOD(I1,3)
         I3 = 1 + MOD(I2,3)
C
C SUBJECT THIS AXES COMBINATION TO ALL NT SYMMETRY OPERATIONS
C
  400    CONTINUE
         DO 420 IT = 1,NT
           IF(IT.EQ.1) GO TO 420
           NN = 9*(IT-1)
           DO 410 I = 1,NAT
              XOLD = CORD(I,I1)
              YOLD = CORD(I,I2)
              ZOLD = CORD(I,I3)
              CALL TRANS(NN)
              CORD2(I,1) = XNEW
              CORD2(I,2) = YNEW
              CORD2(I,3) = ZNEW
  410      CONTINUE
C
C COMPARE GEOMETRY CORD TO CORD2, IF IDENTICAL IANS.EQ.0
C
            CALL EQVATM(AZNUC,CORD,CORD2,IANS,NAT,I1,I2,I3)
C
C IF IANS.NE.0 GEOMETRIES DIDN'T MATCH, AND WE TRY AND
C SWITCH THE LAST TWO AXES.
C
            IF(IANS.NE.0.AND.IICNT.LT.2) THEN
              I3 = 1 + MOD(I1,3)
              I2 = 1 + MOD(I3,3)
              IICNT = IICNT + 1
              GO TO 400
            END IF
C
C THAT DIDN'T WORK LETS TRY A WHOLE NEW COMBINATION
C
            IF(IANS.NE.0.AND.IICNT.EQ.2) GO TO 430
  420    CONTINUE
C
C IF WE GOT HERE IANS.EQ.0 FOR ALL NT SYMMETRY OPERATIONS
C AND WE CONCLUDE THAT THE AXES ARE CORRECT; LETS QUIT.
C
         IF(IANS.EQ.0) GO TO 450
  430 CONTINUE
C
  450 CONTINUE
      DO 460 I = 1,NAT
        CORD2(I,1) = CORD(I,I1)
        CORD2(I,2) = CORD(I,I2)
        CORD2(I,3) = CORD(I,I3)
        CORD(I,1) = CORD2(I,1)
        CORD(I,2) = CORD2(I,2)
        CORD(I,3) = CORD2(I,3)
  460 CONTINUE
C
C     ----- DELETE NON-UNIQUE ATOMS, IF ANY -----
C
      DO 560 I=1,NAT
         XOLD = CORD(I,1)
         YOLD = CORD(I,2)
         ZOLD = CORD(I,3)
         DO 550 IT = 1,NT
            IF(IT.EQ.1) GO TO 550
            NN = 9*(IT-1)
            CALL TRANS(NN)
            DO 510 IAT = I+1,NAT
               TEST = (XNEW-CORD(IAT,1))**2
     *               +(YNEW-CORD(IAT,2))**2
     *               +(ZNEW-CORD(IAT,3))**2
               IF(TEST.LE.TM8) GO TO 520
  510       CONTINUE
C
C      NO NEW ATOM WAS GENERATED, ON TO THE NEXT OPERATION
C
            GO TO 550
C
C      ATOM # IAT WAS GENERATED AND CAN THUS BE DELETED.
C      CONSEQUENTLY NAT DECREASES BY ONE.
C
  520       CONTINUE
            NAT = NAT - 1
            DO 540 J = IAT,NAT
               AATOM(J)  = AATOM(J+1)
               AZNUC(J)  = AZNUC(J+1)
               CORD(J,1) = CORD(J+1,1)
               CORD(J,2) = CORD(J+1,2)
               CORD(J,3) = CORD(J+1,3)
  540       CONTINUE
  550    CONTINUE
  560 CONTINUE
C
C     THE CORD ARRAY NOW HOLDS THE SYMMETRY UNIQUE ATOMS ONLY.
C
      DO 610 I=1,NAT
         UATOM(I) = AATOM(I)
         UZNUC(I) = AZNUC(I)
         XNEW = CORD(I,1)
         YNEW = CORD(I,2)
         ZNEW = CORD(I,3)
         CALL ROT
         CUNIQ(I,1) = XP
         CUNIQ(I,2) = YP
         CUNIQ(I,3) = ZP
  610 CONTINUE
      UATOM(NAT+1) = ENDWRD
      RETURN
C
  810 CONTINUE
      IF (MASWRK) WRITE (IW,9230) MXATM
      CALL ABRT
      STOP
C
 9020 FORMAT(/1X,'THE MOMENTS OF INERTIA ARE (AMU-ANGSTROM**2)'/
     *        1X,'IXX=',F10.3,'   IYY=',F10.3,'   IZZ=',F10.3)
 9230 FORMAT(1X,'*** ERROR, NO MORE THAN',I6,' TOTAL ATOMS')
      END
C*MODULE INPUTC   *DECK NUMERC
      LOGICAL FUNCTION NUMERC(ITEM)
      CHARACTER*8 ITEM
      CHARACTER*1 END
      DIMENSION END(11)
      DATA END/'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '-'/
C
      NUMERC = .FALSE.
      DO 10 I=1,11
        IF(ITEM(1:1).EQ.END(I))  NUMERC = .TRUE.
   10 CONTINUE
      RETURN
      END
C*MODULE INPUTC   *DECK PRAXIS
      SUBROUTINE PRAXIS(C,ZNUC,NAT,NDIM,IXX,IYY,IZZ)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION IXX,IYY,IZZ
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION C(NDIM,3),ZNUC(NDIM)
      DIMENSION AMS(106),ROT(3,3),EIG(3),TRI(6),WRK(3,8),IWRK(3)
C
      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TM8=1.0D-08)
C
C     ----- GENERATE PRINCIPAL AXIS COORDINATES -----
C     THESE ARE RETURNED IN -C-, DESTROYING THE INITIAL SET.
C     ALSO, THE THREE MOMENTS OF INERTIA ARE RETURNED.
C
      CALL RAMS(AMS,2)
C
C     FIRST TRANSLATE TO CENTER OF MASS COORDINATES
C
      TM =ZERO
      XCM=ZERO
      YCM=ZERO
      ZCM=ZERO
      DO 100 I=1,NAT
         NUCZ=INT(ZNUC(I))
         IF(NUCZ.EQ.0) GO TO 100
         TEMP=AMS(NUCZ)
         TM = TM+TEMP
         XCM=XCM+TEMP*C(I,1)
         YCM=YCM+TEMP*C(I,2)
         ZCM=ZCM+TEMP*C(I,3)
  100 CONTINUE
C
      XCM=XCM/TM
      YCM=YCM/TM
      ZCM=ZCM/TM
      IF(ABS(XCM).LT.TM8) XCM=ZERO
      IF(ABS(YCM).LT.TM8) YCM=ZERO
      IF(ABS(ZCM).LT.TM8) ZCM=ZERO
C
      DO 120 I=1,NAT
         C(I,1)=C(I,1)-XCM
         C(I,2)=C(I,2)-YCM
         C(I,3)=C(I,3)-ZCM
  120 CONTINUE
C
C     THEN COMPUTE THE MOMENT OF INERTIA TENSOR.
C
      XX=ZERO
      YY=ZERO
      ZZ=ZERO
      XY=ZERO
      XZ=ZERO
      YZ=ZERO
      DO 200 I=1,NAT
         NUCZ=INT(ZNUC(I))
         IF(NUCZ.EQ.0) GO TO 200
         TEMP=AMS(NUCZ)
         XC=C(I,1)
         YC=C(I,2)
         ZC=C(I,3)
         XX=XX+TEMP*YC*YC+TEMP*ZC*ZC
         YY=YY+TEMP*XC*XC+TEMP*ZC*ZC
         ZZ=ZZ+TEMP*XC*XC+TEMP*YC*YC
         XY=XY-TEMP*XC*YC
         XZ=XZ-TEMP*XC*ZC
         YZ=YZ-TEMP*YC*ZC
  200 CONTINUE
      IF(ABS(XY).GT.TM8) GO TO 220
      IF(ABS(XZ).GT.TM8) GO TO 220
      IF(ABS(YZ).GT.TM8) GO TO 220
      GO TO 500
C
C     DIAGONALIZATION OF THE MOMMENT OF INERTIA TENSOR
C     YIELDS THE ROTATION INTO THE PRINCIPAL AXES.
C
  220 CONTINUE
      TRI(1) = XX
      TRI(2) = XY
      TRI(3) = YY
      TRI(4) = XZ
      TRI(5) = YZ
      TRI(6) = ZZ
C
C        THE JACOBI ROUTINE DOES A BETTER JOB OF GIVING
C        WELL-ORIENTED AXES IN CASE THERE IS A DEGENERATE
C        MOMENT OF INERTIA, SO WE ALWAYS USE THAT ROUTINE.
C
      KDGSAV = KDIAG
      KDIAG = 3
      CALL GLDIAG(3,3,3,TRI,WRK,EIG,ROT,IERR,IWRK)
      KDIAG = KDGSAV
C
C     ADJUST ROTATION MATRIX SO HANDEDNESS IS UNCHANGED
C     A PROPER ROTATION MUST HAVE A POSITIVE DETERMINATE
C
      COF11 = ROT(2,2)*ROT(3,3) - ROT(2,3)*ROT(3,2)
      COF12 = ROT(2,1)*ROT(3,3) - ROT(2,3)*ROT(3,1)
      COF13 = ROT(2,1)*ROT(3,2) - ROT(2,2)*ROT(3,1)
      DET = ROT(1,1)*COF11 - ROT(1,2)*COF12 + ROT(1,3)*COF13
      IF(DET.LT.ZERO) CALL DSCAL(3,-ONE,ROT(1,1),1)
C
C     ROTATE THE COORDINATES TO P.A. FRAME
C
      DO 340 I=1,NAT
         TEMP1=C(I,1)
         TEMP2=C(I,2)
         TEMP3=C(I,3)
         C(I,1)=TEMP1*ROT(1,1) +TEMP2*ROT(2,1) +TEMP3*ROT(3,1)
         C(I,2)=TEMP1*ROT(1,2) +TEMP2*ROT(2,2) +TEMP3*ROT(3,2)
         C(I,3)=TEMP1*ROT(1,3) +TEMP2*ROT(2,3) +TEMP3*ROT(3,3)
  340 CONTINUE
C
      DO 360 J=1,3
         DO 350 I=1,NAT
            IF(ABS(C(I,J)).LT.TM8) C(I,J)=ZERO
  350    CONTINUE
  360 CONTINUE
C
C     TEST FOR DIAGONAL PRODUCT OF INERTIA TENSOR
C
  500 CONTINUE
      XY=ZERO
      XZ=ZERO
      YZ=ZERO
      XX=ZERO
      YY=ZERO
      ZZ=ZERO
      DO 600 I=1,NAT
         NUCZ=INT(ZNUC(I))
         IF(NUCZ.EQ.0) GO TO 600
         TEMP=AMS(NUCZ)
         XC=C(I,1)
         YC=C(I,2)
         ZC=C(I,3)
         XY=XY+TEMP*XC*YC
         XZ=XZ+TEMP*XC*ZC
         YZ=YZ+TEMP*YC*ZC
         XX=XX+TEMP*XC*XC
         YY=YY+TEMP*YC*YC
         ZZ=ZZ+TEMP*ZC*ZC
  600 CONTINUE
C                          A VERY UNLIKELY ERROR CHECK.
      IF(XY.GT.TM8  .OR.  XZ.GT.TM8  .OR.  YZ.GT.TM8) THEN
         IF (MASWRK) WRITE(6,*) 'UNABLE TO GENERATE PRINCIPAL AXES'
         CALL ABRT
         STOP
      END IF
C
C     FINALLY, GENERATE THE MOMENTS OF INERTIA.
C     THE UNITS ARE AMU*DIST**2, WHERE DIST IS THE UNITS -C- HAS.
C
      IXX = YY + ZZ
      IYY = XX + ZZ
      IZZ = XX + YY
      RETURN
      END
C*MODULE INPUTC   *DECK RDMPC
      SUBROUTINE RDMPC(BL,ALPHA,BETA,IZ1,IZ2,IZ3,I)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ZMTALT/ NZMAT2,NZVAR2,NVAR2,NZMTRD,ICOORD
C
      IZ1 = 0
      IZ2 = 0
      IZ3 = 0
      BL = 0.0D+00
      ALPHA = 0.0D+00
      BETA = 0.0D+00
C
      NVAL = 3
      IF(I.EQ.1) NVAL=0
      IF(I.EQ.2) NVAL=1
      IF(I.EQ.3) NVAL=2
      NZVAR2 = NZVAR2+NVAL
C
      IERR = 0
      BL = RFIND('BL      ',IERR)
      IF(IERR.NE.0) CALL ABRT
      IF(I.EQ.2) THEN
         IZ1 = 1
         RETURN
      END IF
      IDUM = IFIND('IDUM1   ',IERR)
      IF(IERR.NE.0) CALL ABRT
C
      ALPHA = RFIND('ALPHA   ',IERR)
      IF(IERR.NE.0) CALL ABRT
      IDUM = IFIND('IDUM2   ',IERR)
      IF(IERR.NE.0) CALL ABRT
C
      BETA = RFIND('BETA    ',IERR)
      IF(IERR.NE.0) CALL ABRT
      IDUM = IFIND('IDUM3   ',IERR)
      IF(IERR.NE.0) CALL ABRT
C
      IZ1 = IFIND('IZ1     ',IERR)
      IF(IERR.NE.0) CALL ABRT
      IZ2 = IFIND('IZ2     ',IERR)
      IF(IERR.NE.0) CALL ABRT
      IZ3 = IFIND('IZ3     ',IERR)
      IF(IERR.NE.0) CALL ABRT
C
      IF(I.EQ.3  .AND.  IZ1+IZ2.EQ.0) THEN
         IZ1=2
         IZ2=1
         IZ3=0
         RETURN
      END IF
C
      IF(I.GT.3  .AND.
     *  (IZ1.EQ.IZ2  .OR.  IZ1.EQ.IZ3  .OR.  IZ2.EQ.IZ3)) THEN
          IF(MASWRK) WRITE(IW,9000) I,IZ1,IZ2,IZ3
          CALL ABRT
          IF(IDUM.NE.0) STOP
      END IF
      RETURN
C
 9000 FORMAT(1X,'ERROR DEFINING ATOM',I4/
     *       1X,'ITS REFERENCE ATOMS',3I4,' CONTAIN A DUPLICATE')
      END
C*MODULE INPUTC   *DECK RDZMT
      SUBROUTINE RDZMT(NATS,IC,AITM,AITM2,NEG)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000)
C
      LOGICAL NUMERC,GOPARR,DSKWRK,MASWRK
      CHARACTER*8 AITM2
      CHARACTER*9 AITM3
      CHARACTER*1 PART1
      CHARACTER*4 TEMPAT
      CHARACTER*4 ATM
C
      EXTERNAL REREAD
C
      COMMON /INFO  / C(MXATM,3),IZAN(MXATM),NATOMS,IUNTRD,ATM(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ZMTALT/ NZMAT2,NZVAR2,NVAR2,NZMTRD,ICOORD
C
      IERR = 0
      IC = 0
      NEG = 1
C
      NZVAR2 = NZVAR2+1
C
      TEMPAT='    '
      LENSTR=-4
      CALL GSTRNG(TEMPAT,LENSTR)
      IF(NUMERC(TEMPAT)) THEN
        CALL REREAD
        IC=IFIND('IZ      ',IERR)
      ELSE
        K=0
        DO 10 I=1,NATS
          IF(TEMPAT.EQ.ATM(I))  K = I
   10   CONTINUE
        IF (K.EQ.0) THEN
           IF (MASWRK) WRITE(IW,9020) TEMPAT,NATS
           CALL ABRT
           STOP
        END IF
        IC = K
      ENDIF
C
C       READ EITHER A NUMBER, OR PERHAPS A STRING
C
      IERR=30
      TEMP=RFIND('AITM    ',IERR)
      IF(IERR.EQ.0) THEN
         AITM = TEMP
      ELSE
         CALL REREAD
         AITM3='         '
         LENSTR=-9
         CALL GSTRNG(AITM3,LENSTR)
         READ(UNIT=AITM3,FMT='(A1,A8)') PART1,AITM2
         IF(PART1.EQ.'-') NEG = -1
         IF(NEG.EQ.-1) RETURN
         READ(UNIT=AITM3,FMT='(A8,A1)') AITM2,PART1
      END IF
      RETURN
C
 9020 FORMAT(1X,'UNKNOWN CONNECTION ATOM=',A4,' ON ZMAT CARD',I4)
      END
C*MODULE INPUTC   *DECK RELVEC
      SUBROUTINE RELVEC(R,J,K)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000)
      CHARACTER*4 ATM
      DIMENSION R(3)
      COMMON /INFO  / C(MXATM,3),IZAN(MXATM),NATOMS,IUNTRD,ATM(MXATM)
C
C     R IS THE RELATIVE VECTOR FROM ATOM K TO ATOM J
C
      DO 10 I=1,3
         R(I)=C(J,I)-C(K,I)
   10 CONTINUE
      RETURN
      END
C*MODULE INPUTC   *DECK RESHUF
      SUBROUTINE RESHUF
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000, TM6=1.0D-06, ONE=1.0D+00)
      PARAMETER (UNITS=ONE/0.52917724924D+00)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      CHARACTER*10 AATOM
C
      COMMON /CONNEC/ IZ(MXATM,4),IATMOR(MXATM),NAT,NRATM,NDATM
      COMMON /INFOA / INAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ZMINF / CORD(MXATM,3),AZNUC(MXATM),CORD2(MXATM,3),
     *                AATOM(MXATM)
C
CJAN
      DO 10 I = 1,NAT
        DO 50 J = 1,3
           CORD2(I,J) = UNITS*CORD2(I,J)
   50   CONTINUE
   10 CONTINUE
C
      DO 200 I = 1,NAT
         DO 100 J = 1,NAT
            TEST = (C(1,I)-CORD2(J,1))**2+(C(2,I)-CORD2(J,2))**2+
     *             (C(3,I)-CORD2(J,3))**2
            IF(TEST.LE.TM6) GO TO 150
  100    CONTINUE
         IF (MASWRK) WRITE(IW,9000) I
  150    CONTINUE
         IATMOR(J) = I
  200 CONTINUE
C
 9000 FORMAT(1X,'GENERATED ATOM ',I4,'NOT FOUND IN ORIGINAL INPUT')
C
      RETURN
      END
C*MODULE INPUTC   *DECK STRIP
      CHARACTER*4 FUNCTION STRIP(ATOM)
      CHARACTER*4 ATOM
      CHARACTER*1 END
      DIMENSION END(11)
      DATA END/'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ' '/
C
      STRIP='    '
      DO 20 I=1,4
        DO 10 J=1,11
          IF(ATOM(I:I).EQ.END(J)) GOTO 30
   10   CONTINUE
        STRIP(I:I)=ATOM(I:I)
   20 CONTINUE
   30 CONTINUE
      RETURN
      END
C*MODULE INPUTC   *DECK UNIVEC
      SUBROUTINE UNIVEC(AL,R)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION AL(3),R(3)
      DATA ONE/1.0D+00/
C
C     AL IS THE UNIT VECTOR OF R
C
      FAC = ONE/SQRT(R(1)*R(1) + R(2)*R(2) + R(3)*R(3))
      DO 10 I=1,3
         AL(I) = R(I) * FAC
   10 CONTINUE
      RETURN
      END
C*MODULE INPUTC   *DECK VECPRD
      SUBROUTINE VECPRD(VP,X,Y)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VP(3),X(3),Y(3)
C
C     VP IS THE VECTOR PRODUCT OF X CROSSED Y
C
      VP(1)=X(2)*Y(3)-X(3)*Y(2)
      VP(2)=X(3)*Y(1)-X(1)*Y(3)
      VP(3)=X(1)*Y(2)-X(2)*Y(1)
      RETURN
      END
C*MODULE INPUTC   *DECK ZMATAB
      SUBROUTINE ZMATAB(MOPAC,NAT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000)
      LOGICAL MOPAC,GOPARR,DSKWRK,MASWRK
      CHARACTER*10 AATOM
      CHARACTER*8  BLANK
      CHARACTER*8  BL2(MXATM), ALPHA2(MXATM), BETA2(MXATM), VAR
      CHARACTER*4  ATMLAB(106), ATM, STRIP
      DIMENSION BL(MXATM),ALPHA(MXATM),BETA(MXATM),A(MXATM),
     *          B(MXATM),D(MXATM),V1(3),V2(3),V3(3),VP(3),VJ(3),
     *          AL1(3),AL2(3),AL3(3),AL4(3),NEG(MXATM,3)
C
      COMMON /CONNEC/ IZ(MXATM,4),IATMOR(MXATM),IQ,NRATM,NDATM
      COMMON /INFO  / C(MXATM,3),IZAN(MXATM),NATOMS,IUNTRD,ATM(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ZMINF / CORD(MXATM,3),AZNUC(MXATM),CORD2(MXATM,3),
     *                AATOM(MXATM)
C
      DATA ANGS/0.52917724924D+00/
      DATA ZERO,ONE,TWO,TM5/0.0D+00,1.0D+00,2.0D+00,1.0D-05/
      DATA PI/3.14159265359D+00/ , ONE80/180.0D+00/
      DATA BADVAL/1234.56578D+09/
      DATA BLANK/'        '/
C
      NERR = 0
      DFACT = ONE
      IF(IUNTRD.EQ.-1) DFACT = ANGS
C
C          * * * * MODEL BUILDER * * * *
C     ORIGINAL AUTHOR - MANY, MANY MOONS AGO:
C                 MARK S. GORDON
C             DEPARTMENT OF CHEMISTRY
C          NORTH DAKOTA STATE UNIVERSITY
C            FARGO, NORTH DAKOTA 58105
C
      DO 110 I=1,MXATM
         BL(I)   = BADVAL
         ALPHA(I)= BADVAL
         BETA(I) = BADVAL
         BL2(I)   = BLANK
         ALPHA2(I)= BLANK
         BETA2(I) = BLANK
         C(I,1)=ZERO
         C(I,2)=ZERO
         C(I,3)=ZERO
  110 CONTINUE
      CALL SETLAB(3,ATMLAB)
C
C        READ THE Z MATRIX INPUT (DUMMY ATOMS HAVE CHARGE 0)
C
      NRATM=0
      NDATM=0
      I=0
      IEOF = 0
  120 CONTINUE
      IF(MOPAC) THEN
         CALL RDCARD('$DATA 8M',IEOF)
      ELSE
         CALL RDCARD('$DATA 8G',IEOF)
      END IF
      IF(IEOF.EQ.1) CALL ABRT
C
      I = I + 1
      ATM(I)='    '
      LENSTR=-4
      CALL GSTRNG(ATM(I),LENSTR)
      IF(ATM(I).EQ.'    ') GO TO 130
      IF(ATM(I).EQ.'$END') THEN
         NATOMS = I - 1
         GO TO 135
      END IF
C
C         READ THE MOPAC OR GAUSS STYLE Z-MATRIX
C
      IF(MOPAC) THEN
         CALL RDMPC(BL(I),ALPHA(I),BETA(I),IZ(I,1),
     *              IZ(I,2),IZ(I,3),I)
      ELSE
         IF(I.GT.1) CALL RDZMT(I,IZ(I,1),BL(I),BL2(I),NEG(I,1))
         IF(I.GT.2) CALL RDZMT(I,IZ(I,2),ALPHA(I),ALPHA2(I),NEG(I,2))
         IF(I.GT.3) THEN
            CALL RDZMT(I,IZ(I,3),BETA(I),BETA2(I),NEG(I,3))
            IERR = 0
            IZ(I,4)=IFIND('IZ',IERR)
            IF (IZ(I,1).EQ.IZ(I,2)  .OR.  IZ(I,1).EQ.IZ(I,3)  .OR.
     *          IZ(I,2).EQ.IZ(I,3)) THEN
               IF(MASWRK) WRITE(IW,9000) I,IZ(I,1),IZ(I,2),IZ(I,3)
               CALL ABRT
            END IF
         END IF
      END IF
C
      BL(I) = DFACT * BL(I)
C
C         TRANSLATE CHEMICAL SYMBOL TO ATOMIC NUMBER
C
      DO 125 J=1,106
         IF(STRIP(ATM(I)).EQ.ATMLAB(J)) GO TO 127
  125 CONTINUE
      IF (MASWRK) WRITE(IW,9050) ATM(I)
      CALL ABRT
      RETURN
C
C          GHOST ATOMS HAVE CHARGE=0
C          NRATM/NDATM = NUMBER OF REAL/DUMMY ATOMS
C
  127 CONTINUE
      IZAN(I) = J
      IF(IZAN(I).GT.104) IZAN(I) = 0
      IF(IZAN(I).GT.0) THEN
         NRATM=NRATM+1
      ELSE
         NDATM=NDATM+1
      END IF
      GO TO 120
C
C        DONE PROCESSING Z-MATRIX SECTION
C
  130 CONTINUE
      NATOMS = I-1
      IF(MOPAC) GO TO 135
C
C        READ VALUE FOR ANY SYMBOLS WHICH MAY BE IN GAUSSIAN Z-MATRIX
C
  131 CONTINUE
      CALL RDCARD('$DAT 10G',IEOF)
      IF(IEOF.EQ.1) CALL ABRT
      VAR='        '
      LENSTR=-8
      CALL GSTRNG(VAR,LENSTR)
      IF(VAR.EQ.BLANK) GO TO 131
      IF(VAR.EQ.'$END    ') GO TO 133
      AVAR = RFIND('AVAR   ',IERR)
      IF(IERR.EQ.1) CALL ABRT
C
      MATCH = 0
      DO 132 I=1,NATOMS
         IF(   BL2(I).EQ.VAR) THEN
            BL(I) = NEG(I,1)*AVAR
            MATCH = MATCH+1
         END IF
         IF(ALPHA2(I).EQ.VAR) THEN
            ALPHA(I) = NEG(I,2)*AVAR
            MATCH = MATCH+1
         END IF
         IF( BETA2(I).EQ.VAR) THEN
            BETA(I) = NEG(I,3)*AVAR
            MATCH = MATCH+1
         END IF
  132 CONTINUE
      IF(MATCH.GT.0) GO TO 131
C
C         ONE OF THE VARIABLE STRINGS WAS NOT GIVEN IN THE Z-MATRIX
C
      IF (MASWRK) WRITE(IW,9030) VAR
      NERR=NERR+1
      GO TO 131
C
C         MAKE SURE ALL SYMBOLS RECEIVED A VALUE
C
  133 CONTINUE
      DO 134 I=1,NATOMS
         IF(   BL2(I).NE.BLANK .AND.    BL(I).EQ.BADVAL) THEN
            IF (MASWRK) WRITE(IW,9035) BL2(I)
            NERR=NERR+1
         END IF
         IF(ALPHA2(I).NE.BLANK .AND. ALPHA(I).EQ.BADVAL) THEN
            IF (MASWRK) WRITE(IW,9035) ALPHA2(I)
            NERR=NERR+1
         END IF
         IF( BETA2(I).NE.BLANK .AND.  BETA(I).EQ.BADVAL) THEN
            IF (MASWRK) WRITE(IW,9035) BETA2(I)
            NERR=NERR+1
         END IF
  134 CONTINUE
C
C        PRINT OUT NUMERIC Z MATRIX
C
  135 CONTINUE
      IF(NERR.GT.0) THEN
         IF (MASWRK) WRITE(IW,*)
     *      'RESUBMIT AFTER FIXING YOUR Z-MATRIX ERROR(S)'
         CALL ABRT
         STOP
      END IF
      IF (MASWRK) THEN
      WRITE(IW,9040)
      DO 137 I=1,NATOMS
         IF(I.EQ.1) WRITE(IW,9055) ATM(I)
         IF(I.EQ.2) WRITE(IW,9055) ATM(I),IZ(I,1),BL(I)
         IF(I.EQ.3) WRITE(IW,9055) ATM(I),IZ(I,1),BL(I),
     *                                    IZ(I,2),ALPHA(I)
         IF(I.GE.4) WRITE(IW,9055) ATM(I),IZ(I,1),BL(I),
     *                                    IZ(I,2),ALPHA(I),
     *                                    IZ(I,3),BETA(I),IZ(I,4)
  137 CONTINUE
      ENDIF
C
      M=NATOMS
      IF (M.LT.3) GO TO 180
      DO 140 J=3,M
         ALPHA(J) = ALPHA(J)*PI/ONE80
  140 CONTINUE
      IF(M.LT.4) GO TO 180
      DO 160 J=4,M
         BETA(J) = BETA(J)*PI/ONE80
  160 CONTINUE
  180 CONTINUE
C
C              HANDLE ATOMS 1, 2, 3 AS SPECIAL CASES
C
      C(1,1)=ZERO
      C(1,2)=ZERO
      C(1,3)=ZERO
C
      C(2,1)=ZERO
      C(2,2)=ZERO
      C(2,3)=BL(2)
      IF(M.LT.3) GO TO 700
C
      C(3,1)=BL(3)*SIN(ALPHA(3))
      C(3,2)=ZERO
      IF(IZ(3,1).EQ.1) C(3,3) = BL(3)*COS(ALPHA(3))
      IF(IZ(3,1).EQ.2) C(3,3) = C(2,3)-BL(3)*COS(ALPHA(3))
      I=3
      IF(ABS(C(3,1)).GE.TM5) GO TO 240
      DO 220 I = 4,MXATM
         IF(I.LE.M.AND.ABS(C(I-1,1)).LT.TM5) GO TO 210
            GO TO 240
  210    C(I,1) = BL(I)*SIN(ALPHA(I))
         C(I,2) = ZERO
         ITEMP = IZ(I,1)
         JTEMP = IZ(I,2)
         C(I,3)=C(ITEMP,3)-BL(I) * COS(ALPHA(I)) *
     *           SIGN(ONE,C(ITEMP,3)-C(JTEMP,3))
  220 CONTINUE
  240 CONTINUE
      K=I
      IF(I.EQ.3) K=4
C
C               LOOP OVER COORDINATES OF ATOMS 4 TO M
C
      IF(K.GT.M) GO TO 700
      DO 600 J = K,M
      IF(IZ(J,4).NE.0) GO TO 400
C
C                  ****** USE DIHEDRAL ANGLE *****
C
      CALL RELVEC(V1,IZ(J,2),IZ(J,3))
      CALL UNIVEC(AL1,V1)
C
C     V1 IS A VECTOR FROM ATOM IZ(J,3) TO ATOM IZ(J,2)
C     AL1 IS A UNIT VECTOR OF V1.
C     BL(IZ(J,2)) IS BONDLENGTH OF ATOM IZ(J,2)
C
      CALL RELVEC(V2,IZ(J,1),IZ(J,2))
      CALL UNIVEC(AL2,V2)
C
C     V2 IS THE VECTOR FROM ATOM IZ(J,2) TO ATOM IZ(J,1)
C     AL2 IS THE UNIT VECTOR OF V2.
C     BL(IZ(J,1)) IS BONDLENGTH DEFINING ATOM IZ(J,1)
C
      CALL VECPRD(VP,AL1,AL2)
      DO 300 I = 1,3
         AL3(I) = VP(I)/
     *        SQRT(ONE-(AL1(1)*AL2(1)+AL1(2)*AL2(2)+AL1(3)*AL2(3))**2)
  300 CONTINUE
C
C     AL3 IS THE UNIT VECTOR OF VP.ALPHA(IZ(J,1))IS THE ANGLE DEFINING
C     ATOM IZ(J,1)
C
      CALL VECPRD(AL4,AL3,AL2)
C
C     AL4,AL3,AL2 ARE A NEW SET OF MUTUALLY ORTHOGONAL AXES. WE WILL NOW
C     GET COORD OF ATOM J IN TERMS OF THESE AXES RELATIVE TO ATOM IZ(J)
C
      DO 320 I=1,3
         VJ(I) = BL(J)*(-AL2(I)*COS(ALPHA(J))+AL4(I)*SIN(ALPHA(J))*
     *           COS(BETA(J))+AL3(I)*SIN(ALPHA(J))*SIN(BETA(J)))
C
C     VJ IS THE VECTOR FROM ATOM IZ(J,1) TO ATOM J
C
         ITEMP = IZ(J,1)
         C(J,I) = VJ(I) + C(ITEMP,I)
         IF(ABS(C(J,I)).LT.TM5) C(J,I)=ZERO
  320 CONTINUE
      GO TO 600
C
C              ***** USE SECOND BOND ANGLE *****
C
  400 CONTINUE
      IF(IABS(IZ(J,4)).NE.1) GO TO 500
      CALL RELVEC(V1,IZ(J,1),IZ(J,3))
      CALL UNIVEC(AL1,V1)
      CALL RELVEC(V2,IZ(J,2),IZ(J,1))
      CALL UNIVEC(AL2,V2)
      ZETA = -(AL1(1)*AL2(1)+AL1(2)*AL2(2)+AL1(3)*AL2(3))
      OMZSQ = ONE - ZETA * ZETA
      A(J) = (-COS(BETA(J))+ZETA*COS(ALPHA(J)))/OMZSQ
      B(J) = (COS(ALPHA(J))-ZETA*COS(BETA(J)))/OMZSQ
      TEMP=ZERO
      IF(ZETA.LT.ZERO) TEMP=PI
      GAMMA = PI/TWO
      IF(ZETA.NE.ZERO) GAMMA = ATAN (SQRT(OMZSQ)/ZETA) +TEMP
      IF (ABS(GAMMA+ALPHA(J)+BETA(J)-TWO*PI).GT.TM5)
     *     D(J) = IZ(J,4)*(SQRT(ONE+A(J)*COS(BETA(J))-B(J)*
     *            COS(ALPHA(J))))/SQRT(OMZSQ)
      CALL VECPRD(V3,AL1,AL2)
      DO 420 I = 1,3
         AL3(I) = A(J)*AL1(I)+B(J)*AL2(I)+D(J)*V3(I)
         VJ(I) = BL(J)*AL3(I)
         ITEMP = IZ(J,1)
         C(J,I) = VJ(I)+C(ITEMP,I)
         IF(ABS(C(J,I)).LT.TM5) C(J,I)=ZERO
  420 CONTINUE
      GO TO 600
C
C            ***** USE BOND-PLANE ANGLE ******
C
C     THIS PART IS USED IF ATOM IS DEFINED BY ONE BONDANGLE AND BY THE
C     ANGLE WHICH BOND IZ(J,1)-J MAKES WITH THE PLANE OF IZ(J,1),IZ(J,2)
C     AND IZ(J,3)
C
  500 CONTINUE
      CALL RELVEC(V1,IZ(J,1),IZ(J,3))
      CALL UNIVEC(AL1,V1)
      CALL RELVEC(V2,IZ(J,2),IZ(J,1))
      CALL UNIVEC(AL2,V2)
      ZETA = -(AL1(1)*AL2(1)+AL1(2)*AL2(2)+AL1(3)*AL2(3))
      OMZSQ = ONE - ZETA * ZETA
      CALL VECPRD(V3,AL1,AL2)
C
C     BETA(J) HERE REFERS TO 90-ANGLE WHICH BOND IZ(J,1)-J MAKES
C     WITH THE AFOREMENTIONED PLANE
C
      V3MAG = SQRT(V3(1)*V3(1) + V3(2)*V3(2) + V3(3)*V3(3))
      A(J)=(COS(BETA(J))/OMZSQ)*V3MAG
      B(J)=SQRT((ONE-COS(ALPHA(J))**2
     *              - A(J)*COS(BETA(J))*V3MAG)/OMZSQ)
      IF(IZ(J,4).EQ.-2) B(J) = -B(J)
      D(J) = B(J)*ZETA+COS(ALPHA(J))
      DO 520 I=1,3
         AL3(I) = B(J)*AL1(I)+D(J)*AL2(I)+A(J)*V3(I)
         VJ(I) = BL(J)*AL3(I)
         ITEMP = IZ(J,1)
         C(J,I) = VJ(I) + C(ITEMP,I)
         IF(ABS(C(J,I)).LT.TM5) C(J,I)=ZERO
  520 CONTINUE
  600 CONTINUE
C
  700 CONTINUE
      IQ = 0
      DO 720 J=1,M
         IQ = IQ + 1
         AZNUC(IQ) = IZAN(IQ)
         AATOM(IQ) = ATM(IQ)
         CORD(IQ,1) = C(IQ,1)
         CORD(IQ,2) = C(IQ,2)
         CORD(IQ,3) = C(IQ,3)
  720 CONTINUE
      NAT = M
      RETURN
C
 9000 FORMAT(1X,'ERROR DEFINING ATOM',I4/
     *       1X,'ITS REFERENCE ATOMS',3I4,' CONTAIN A DUPLICATE')
 9030 FORMAT(1X,'A VALUE WAS GIVEN FOR A VARIABLE ',A8,
     *    ' NOT FOUND IN THE Z-MATRIX.')
 9035 FORMAT(1X,'THE SYMBOL ',A8,' GIVEN IN THE Z-MATRIX',
     *    ' DID NOT RECEIVE A VALUE.')
 9040 FORMAT(/1X,'YOUR FULLY SUBSTITUTED Z-MATRIX IS')
 9050 FORMAT(/1X,'ILLEGAL CHEMICAL SYMBOL, ',A4,' USED')
 9055 FORMAT(1X,A4,1X,I3,F12.7,I3,F10.4,I3,F10.4,I3)
      END
