C 21 APR 08 - JAB - GTBOND: FIX FOR SEVERAL TYPES OF FRAGMENTS
C  4 MAR 08 - BN  - GINVR: SKIP BOMBOUT DURING VSCF RUNS
C  7 DEC 07 - JAB - MOVE RRD CALL TO FIX MWS' BUG IN 25 JUN 07 CHANGE
C 25 JUN 07 - CHC - ALLOW FROZEN COORDS IFZMAT/VALUE WHEN AUTO=.FALSE.
C 21 MAY 07 - JMM - GTBOND: SKIP DUMMIES FOR MAKEFP RUNS
C 24 MAR 07 - MWS - PAD FRGINF COMMON BLOCK
C 19 SEP 05 - IA  - SYNCHRONIZE FRGINF COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 14 MAR 05 - HL  - CHANGE FMASS IN FRGMSS COMMON
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK FRGINF
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  4 NOV 03 - JMR - ZMATIN: RAISE ANGLE/TORSION ESTIMATES
C  3 SEP 03 - CHC - ZMAT2: SET CORRECT NAME IO SIZES
C 20 JUN 02 - JAB - INCREASE NUMBER OF LINEAR BENDS PERMITTED
C 22 MAY 02 - MWS - ZMAT2: AVOID VARIABLE COUNT CHECK IF NZVAR=0 INPUT
C 17 APR 02 - MWS - SYNCH UP EFMULT AND FRGINF COMMON
C 24 JAN 02 - KRG - PREVENT USE OF DLC FOR LINEAR SYSTEMS
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAME IO CALL
C 13 JUN 01 - MWS - ZMATIN: RAISE TORSION STORAGE GUESS
C 16 FEB 00 - CHC,MWS - ZMAT2: CHANGE MEMORY ALLOCATIONS
C 29 AUG 99 - CHC - ZMAT2: GEO. CONSTRAINTS WITH QMMM
C 27 SEP 98 - CHC - CONSTRAINT GEOMETRY OPTIMIZATION IN DLC
C 12 SEP 98 - JAB - ZMAT2: MOVE ASSIGNMENT OF -NROW-
C 31 MAR 98 - JRS   PLUGS ADDED FOR AUTOMATIC DLC GENERATION
C 13 MAR 98 - MWS - PZANDG: EXIT W/O PRINTING IF 1,2,3 ATOMS ONLY,
C                   ZMTSYM: LEAVE ALL SYM COORDS TYPE=-1 IF DLC USED
C  6 JAN 98 - JAB - ZMAT2:  SET DECOMP=.F. IF SYM. COORDS ARE USED
C                   SYM. COORD TYPES= -1 (STRETCHES) AND -2 (OTHERS)
C 13 NOV 96 - JRS - ZMATIN,TORS: READ BIG SIJ MATRIX, CHANGE TORS TEST
C 17 OCT 96 - MWS - GTBOND: CORRECT COMMON DECLARATION
C 18 SEP 96 - MWS - GTBOND: ADD CODE FOR FRAGMENT BONDS
C 18 APR 96 - RZ  - ZMATIN: INCREASE DIMENSION FOR IJS/SIJ FROM 7 TO 30
C 26 JUL 95 - PND - ZMATIN: INCREASE DIMENSION FOR IJS/SIJ TO 7
C 19 JUN 95 - MWS - GINVR: DOUBLE DEFINITION OF "SINGULAR".
C  1 FEB 95 - MWS - LIBE: NEED TO SAVE APTS INPUT
C 12 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 17 DEC 93 - MWS - GTBOND: COVALENT RADII DEFINE BONDED ATOM LIST
C 15 DEC 93 - MWS - BANDBI,BMAT,GINVR,GTBOND,ZMAT2: CHANGE DAF NOS.
C                   BMAT,BANDBI,GINVR: SAVE B MATRIX OF REDUNDANT COORDS
C                   ZMTSYM: FIX IZMAT2 GENERATION FOR LINEAR BENDS
C 11 NOV 93 - JHJ - IMPLEMENT USE OF SYMMETRY COORDINATES
C 16 JUL 93 - MWS - GINVR: CHANGE DEBUG PRINTING
C 29 MAR 93 - TLW - MAKE B AND B INVERSE RECTANGULAR
C 18 MAR 93 - JHJ - ZMAT2: PUNCH A TEMPLATE $ZMAT GROUP
C 22 OCT 92 - JAB - FIX SIGN ERRORS IN BDIHA5
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C  3 FEB 92 - TLW - BDIHA5: PUT CALL ABRT AND STOP AT END OF PROCEDURE
C 26 JAN 92 - MWS - FILLZM: CAN'T MAKE $ZMAT IF DUMMIES USED IN $DATA
C  8 JAN 92 - TLW - MAKE WRITES PARALLEL; ADD COMMON PAR
C  7 NOV 91 - MWS - BINVR: PRINT ERROR MESSAGE
C 21 AUG 91 - MWS - PZANDG,ZMAT2: ADD ICOORD TO /ZMTALT/
C 19 AUG 91 - MWS - FILLZM: FIX TEST FOR 2ND BOND ANGLE
C  9 NOV 90 - JHJ - ZMATIN,PZANDG,PRZDAT:ADD FULLY SUBST Z-MAT PRINT,
C                   ADD ROUTINE FILLZM, WHICH IS CALLED BY ZMATIN.
C  3 NOV 90 - MWS - BMAT,GINVR: ALLOW FOR ZERO MASS ATOMS
C  7 MAR 90 - MWS - ADD NEW ROUTINE GTBOND
C 22 SEP 89 - MWS - CHANGE ANGSTROMS-BOHR CONVERSION FACTOR
C 25 FEB 89 - STE - WRITE I80 INSTEAD OF 180 IN BANDBI, TORS PARAMETER
C 30 JAN 89 - MWS - PASS MASSES IN /MASSES/ INSTEAD OF RAMS CALL
C 10 JAN 89 - MR  - ADD SIGN TO DIHEDRAL ANGLES IN TORS
C 22 AUG 88 - JAB - CALL ABRT IN GINVR IF TOO MANY ZERO EIGENVALUES
C 08 AUG 88 - JAB - ADD ROUTINE GINVR FOR VIBRATIONAL ANALYSIS
C 27 JUL 88 - JAB - MAKE NVAR .GE. 3N-6 OK IN ZMATIN IF RUNTYP=FORCE
C 13 JUN 88 - JAB - ADD GTORS ROUTINE FOR "GHOST ATOM" TORSION
C 14 NOV 87 - STE - USE EXETYP; USE PARAMETERS
C 11 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
C 14 NOV 85 - LAM - IMPROVE ERROR MESSAGES IN OPLA
C  1 NOV 85 - STE - ZMATIN: FIX DYNAMIC MEMORY FOR LOADFM.NE.0
C 23 OCT 85 - STE - USE GENERIC ABS,MAX,SIGN,SIN,SQRT
C                   BMAT: MAKE SURE B MATRIX IS COMPLETELY ZEROED
C 11 SEP 85 - STE - ZMATIN: INSERT NVAR=0; DELETE 9028 FORMAT
C                   LIBE:  LIMIT ACOS VALUE RANGE
C                   INITIALIZE MISC VARIABLES, IMPROVE DIAGNOSTICS
C 17 APR 85 - MWS - INTERNAL COORD PRINTOUT IN COL 1-70
C  9 MAR 85 - MWS - FIX LINEAR BENDS, ALLOW 3N-5 INTERNAL COORDS,
C                   PRINT CODED ZMATRIX INPUT
C 24 APR 84 - STE - FIX ANGS/BOHR CONSTANT TYPO IN PRZDAT
C  9 APR 84 - STE - REPLACE BINVR, USE RAMS, DELETE AMASS,MATINV
C  8 MAR 84 - STE - GET CORRECT VALUE OF PI FROM ACOS
C 17 NOV 83 - STE - FIX/DELETE FORMATS IN BANDBI,BINVR,BMAT
C  4 NOV 83 - STE - COMBINE ZGMATP AND ZMATPR INTO PZANDG, CLEAN BMAT
C  6 MAY 83 - MWS - FIX BUG IN PRZDAT
C  7 MAR 83 - MWS - SPLIT AWAY FROM INPUT, WRITE BOHRS AND
C                   RADIANS ON DAF FOR SADPOINT RUNS.
C 25 FEB 83 - STE - ADD  BDIHA5 AND ASSOCIATED CHANGES
C                   IN ZMATIN, BMAT, AND PRZDAT.
C
C*MODULE ZMATRX  *DECK BANDBI
      SUBROUTINE BANDBI
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL LINEAR
C
      PARAMETER (MXATM=2000)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,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 /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
C
C     ----- CALCULATE B AND B INVERSE MATRICES -----
C     AT THIS POINT IT IS GOOD TO REALIZE THAT
C     NVAR  = NUMBER OF COORDS IN THE GEOMETRY SEARCH, 3N-6 OR 3N-5
C     NZVAR = NUMBER OF INPUT INTERNAL COORDINATES, ADJUSTED UPWARDS
C             BY 1 FOR EACH LINEAR BEND.  THUS IT MAY NO LONGER BE
C             EQUAL TO THE USER'S INPUT.  WE HAVE NZVAR.GE.NVAR
C     NRWO  = ??
C
      NCOORD = 3*NAT
C
C     ----- GROW FAST MEMORY -----
C
      CALL VALFM(LOADFM)
      LIZMAT = 1      + LOADFM
      LZMAT  = LIZMAT + NZMAT
      LB     = LZMAT  + NZVAR
      LBINV  = LB     + NZVAR*NCOORD
      LEVEC  = LBINV  + NVAR*NCOORD
      LG     = LEVEC  + NVAR*NVAR
      LRECPM = LG     + NVAR*NVAR
      LEVAL  = LRECPM + NCOORD
      LSCR   = LEVAL  + NVAR
      LIA    = LSCR   + 8*NVAR
      LGLIN  = LIA    + NVAR
      LSYM   = LGLIN  + (NVAR*NVAR+NVAR)/2
      LT     = LSYM   + NZVAR*NZVAR
      LAST   = LT     + NCOORD
      NEED = LAST - LOADFM
      CALL GETFM(NEED)
C
C     ----- GET IZMAT AND SYMMETRY MATRIX FROM DIRECT ACCESS FILE -----
C
      CALL DAREAD(IDAF,IODA,X(LIZMAT),NZMAT,40,1)
      CALL DAREAD(IDAF,IODA,X(LSYM),NZVAR*NZVAR,46,0)
C
C     ----- CALCULATE THE B MATRIX -----
C
      CALL BMAT(X(LB),X(LIZMAT),X(LZMAT),
     *          X(LSYM),X(LT),NZVAR,NVAR,NCOORD)
C
C     ----- CALCULATE THE B INVERSE MATRIX -----
C
      CALL GINVR(X(LB),X(LBINV),X(LEVEC),X(LG),X(LRECPM),X(LEVAL),
     *           X(LSCR),X(LIA),X(LGLIN),NVAR,NCOORD,NAT)
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE ZMATRX  *DECK BDIHA5
C*
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY - USDOE)  DATE: 25 FEB 1983
C*
C*    PURPOSE -
C*       EVALUATE B MATRIX AND DIHEDRAL ANGLE BETWEEN 5 ATOMS FORMING
C*                *            ***      *             *
C*       TWO PLANES HAVING ONE ATOM IN COMMON
C*
C*    ON ENTERING -
C*       IW     - INTEGER (LOGICAL UNIT NO.)
C*                ERROR MESSAGE FILE
C*       NOINT  - INTEGER
C*                NO. OF INTERNAL COORDINATE BEING EVALUATED
C*       NDIM   - INTEGER
C*                NO. OF COORDINATES = 3 * NATOMS
C*       IZMAT  - INTEGER (5)
C*                LIST OF ATOMS INVOLVED
C*       CARTC  - W.P. REAL (3,NATOMS)
C*                CARTESIAN COORDINATES FOR ATOMS
C*
C*    ON EXITING -
C*       DIHANG - W.P. REAL
C*                DIHEDRAL ANGLE IN DEGREES
C*       BMAT   - W.P. REAL (NDIM,3*NAT)
C*                B MATRIX (TRANSFORMS INTERNAL COORD TO CART. COORD
C*                          TO FIRST ORDER)
C*
      SUBROUTINE BDIHA5(IW,NOINT,NDIM,IZMAT,CARTC,DIHANG,BMAT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION IZMAT(5), CARTC(3,1), BMAT(NDIM,1),
     *          A(3),B(3),C(3),D(3),E1(3),E2(3),E3(3)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
C           THE FIRST THREE ATOMS ARE ON THE FIRST PLANE
C           THE  LAST THREE ATOMS ARE ON THE SECOND PLANE
C
      IATOM = IZMAT(1)
      JATOM = IZMAT(2)
      KATOM = IZMAT(3)
      LATOM = IZMAT(4)
      MATOM = IZMAT(5)
C
C           DEFINE VECTORS FROM PIVOT (THIRD) ATOM THAT DEFINE PLANES
C                    A = I - K   AND  B = J - K   DEFINE FIRST PLANE
C                    C = L - K   AND  D = M - K   DEFINE SECOND PLANE
C
      DO 110 IXYZ=1,3
         A(IXYZ) = CARTC(IXYZ,IATOM) - CARTC(IXYZ,KATOM)
         B(IXYZ) = CARTC(IXYZ,JATOM) - CARTC(IXYZ,KATOM)
         C(IXYZ) = CARTC(IXYZ,LATOM) - CARTC(IXYZ,KATOM)
         D(IXYZ) = CARTC(IXYZ,MATOM) - CARTC(IXYZ,KATOM)
  110 CONTINUE
C
C              E1 = A X B  IS PERPENDICULAR TO FIRST PLANE
C
      E1(1) = A(2)*B(3) - A(3)*B(2)
      E1(2) = A(3)*B(1) - A(1)*B(3)
      E1(3) = A(1)*B(2) - A(2)*B(1)
      AXBMAG = SQRT(E1(1)**2 + E1(2)**2 + E1(3)**2)
      IF(ONE+AXBMAG.EQ.ONE) GO TO 810
      AXBMI = ONE / AXBMAG
C
C              E2 = C X D  IS PERPENDICULAR TO SECOND PLANE
C
      E2(1) = C(2)*D(3) - C(3)*D(2)
      E2(2) = C(3)*D(1) - C(1)*D(3)
      E2(3) = C(1)*D(2) - C(2)*D(1)
      CXDMAG = SQRT(E2(1)**2 + E2(2)**2 + E2(3)**2)
      IF(ONE+CXDMAG.EQ.ONE) GO TO 820
      CXDMI = ONE / CXDMAG
C
C              NORMALIZE E1 AND E2
C
      DO 120 IXYZ=1,3
         E1(IXYZ) = AXBMI * E1(IXYZ)
         E2(IXYZ) = CXDMI * E2(IXYZ)
  120 CONTINUE
C
C        THE DIRECTION OF E2 IS DEFINED SUCH THAT WHEN THE PLANES ARE
C        ROTATED ALONG THE INTERSECTION E3 (= E1 X E2) SO THEY COINCIDE
C        (DIHEDRAL ANGLE = 0), THEN ATOMS 1 AND 4 LIE ON THE SAME
C        SIDE OF E3.  IF ATOMS 1 AND/OR 4 LIE ON E3, THEN
C        ATOMS 2 AND/OR 5 ARE USED.
C
      ADOTE2 = A(1)*E2(1) + A(2)*E2(2) + A(3)*E2(3)
      IF(ONE+ADOTE2.EQ.ONE) ADOTE2 = B(1)*E2(1)+B(2)*E2(2)+B(3)*E2(3)
      CDOTE1 = C(1)*E1(1) + C(2)*E1(2) + C(3)*E1(3)
      IF(ONE+CDOTE1.EQ.ONE) CDOTE1 = D(1)*E1(1)+D(2)*E1(2)+D(3)*E1(3)
      IF(ADOTE2*CDOTE1 .LT. ZERO) GO TO 140
C
C        CHANGE SENSE OF E2 SO A (OR B) WILL BE ON SAME SIDE
C        OF E1 X E2 AS C (OR D)
C
         ADOTE2=-ADOTE2
         DO 130 IXYZ=1,3
            E2(IXYZ) = -E2(IXYZ)
  130    CONTINUE
  140 CONTINUE
C
C           EVALUATE DIHEDRAL ANGLE    0 < ABS(DIHANG) < 180
C           ANGLE IS POSITIVE FOR SECOND PLANE "ABOVE" ATOM 1
C           ANGLE IS NEGATIVE FOR SECOND PLANE "BELOW" ATOM 1
C
      E1DE2 = E1(1)*E2(1) + E1(2)*E2(2) + E1(3)*E2(3)
      IF(ABS(E1DE2).GT.ONE) GO TO 840
      IF(ABS(E1DE2).EQ.ONE) GO TO 830
      DIHANG = ACOS(E1DE2)
      IF(ADOTE2 .GT. ZERO) DIHANG = -DIHANG
      SINDI = ONE / SIN(DIHANG)
C
C         E3 = E1 X E2  IS ALONG THE INTERSECTION LINE OF THE TWO PLANES
C
      E3(1) = SINDI * (E1(2)*E2(3) - E1(3)*E2(2))
      E3(2) = SINDI * (E1(3)*E2(1) - E1(1)*E2(3))
      E3(3) = SINDI * (E1(1)*E2(2) - E1(2)*E2(1))
C
      F1 =  AXBMI * (B(1)*E3(1) + B(2)*E3(2) + B(3)*E3(3))
      F2 = -AXBMI * (A(1)*E3(1) + A(2)*E3(2) + A(3)*E3(3))
      F4 =  CXDMI * (D(1)*E3(1) + D(2)*E3(2) + D(3)*E3(3))
      F5 = -CXDMI * (C(1)*E3(1) + C(2)*E3(2) + C(3)*E3(3))
C
      I = 3*IATOM - 3
      J = 3*JATOM - 3
      K = 3*KATOM - 3
      L = 3*LATOM - 3
      M = 3*MATOM - 3
C
C           CALCULATE B MATRIX ELEMENTS
C
      DO 160 N=1,3
         B1 = E1(N) * F1
         B2 = E1(N) * F2
         B4 = E2(N) * F4
         B5 = E2(N) * F5
         B3 = -(B1+B2+B4+B5)
         BMAT(NOINT,I+N) = B1
         BMAT(NOINT,J+N) = B2
         BMAT(NOINT,K+N) = B3
         BMAT(NOINT,L+N) = B4
         BMAT(NOINT,M+N) = B5
  160 CONTINUE
      RETURN
C
C               ERROR MESSAGE SECTION
C
  810 CONTINUE
      IF (MASWRK) WRITE(IW,910)
      GO TO 880
C
  820 CONTINUE
      IF (MASWRK) WRITE(IW,920)
      GO TO 880
C
  830 CONTINUE
      IF (MASWRK) WRITE(IW,930)
      GO TO 880
C
  840 CONTINUE
      IF (MASWRK) WRITE(IW,940)
C
  880 CONTINUE
      IF (MASWRK) THEN
      WRITE(IW,901) NOINT,IATOM,(CARTC(IXYZ,IATOM),IXYZ=1,3)
     *                   ,JATOM,(CARTC(IXYZ,JATOM),IXYZ=1,3)
     *                   ,KATOM,(CARTC(IXYZ,KATOM),IXYZ=1,3)
     *                   ,LATOM,(CARTC(IXYZ,LATOM),IXYZ=1,3)
     *                   ,MATOM,(CARTC(IXYZ,MATOM),IXYZ=1,3)
      WRITE(IW,900)
      END IF
      CALL ABRT
      STOP
C
  900 FORMAT('0THESE MESSAGES BROUGHT TO YOU BY ROUTINE BDIHA5 *****')
  901 FORMAT('0INTERNAL COORDINATE',I4,' (DIHEDRAL ANGLE BETWEEN TWO'
     *      ,' PLANES WITH ONE ATOM IN COMMON) USES THESE ATOMS'/
     *       (I10,3F20.6))
  910 FORMAT('0***** ERROR - FIRST THREE ATOMS ARE LINEAR *****')
  920 FORMAT('0***** ERROR -- LAST THREE ATOMS ARE LINEAR *****')
  930 FORMAT('0***** ERROR - ALL FIVE ATOMS LIE IN THE SAME PLANE ****')
  940 FORMAT('0***** ERROR - COSIN OF DIHEDRAL ANGLE IS GREATER THAN 1')
      END
C*MODULE ZMATRX  *DECK BEND
      SUBROUTINE BEND(EQVAL,NOINT,I,J,K,C,B,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
C     -----THIS ROUTINE COMPUTES THE B MATRIX ELEMENTS OF A
C          VALENCE ANGLE BENDING COORDINATE AS DEFINED BY WILSON.
C          SEE PP. 56-58 OF "MOLECULAR VIBRATIONS" BY E.B. WILSON,
C          J.C. DECIUS AND P.C. CROSS, MCGRAW-HILL, NY, 1955 ----
C
C     -----I AND K ARE THE NUMBERS OF THE END ATOMS.  J IS THE
C          NUMBER OF THE CENTRAL ATOM -----
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION C(3,1),B(NDIM,1)
      DIMENSION RJI(3),RJK(3),EJI(3),EJK(3)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TOL=1.0D-07)
C
C
      PI = ACOS(-ONE)
      DJISQ = ZERO
      DJKSQ = ZERO
C
      DO 100 M = 1,3
         RJI(M) = C(M,I)-C(M,J)
         RJK(M) = C(M,K)-C(M,J)
         DJISQ = DJISQ+RJI(M)*RJI(M)
         DJKSQ = DJKSQ+RJK(M)*RJK(M)
  100 CONTINUE
      DJI = SQRT(DJISQ)
      DJK = SQRT(DJKSQ)
      DOTJ = ZERO
      DO 120 M = 1,3
         EJI(M) = RJI(M)/DJI
         EJK(M) = RJK(M)/DJK
         DOTJ = DOTJ+EJI(M)*EJK(M)
  120 CONTINUE
      IF (ONE.LE.ABS(DOTJ)) GO TO 200
      SINJ = SQRT(ONE-DOTJ*DOTJ)
      NOCOL1 = 3*(I-1)
      NOCOL2 = 3*(K-1)
      NOCOL3 = 3*(J-1)
      DO 160 M = 1,3
         B1 = ((DOTJ*EJI(M)-EJK(M)))/(DJI*SINJ)
         B2 = ((DOTJ*EJK(M)-EJI(M)))/(DJK*SINJ)
         B(NOINT,NOCOL1+M) = B1
         B(NOINT,NOCOL2+M) = B2
         B(NOINT,NOCOL3+M) = -(B1+B2)
  160 CONTINUE
      DOT = ZERO
      DO 180 M = 1,3
  180 DOT = DOT+EJI(M)*EJK(M)
      EQVAL =   ACOS(DOT)
      IF(ABS(EQVAL).LT.TOL) EQVAL=ZERO
      IF(ABS(EQVAL-PI).LT.TOL) EQVAL=PI
      RETURN
C
  200 CONTINUE
      IF (MASWRK) WRITE (IW,9008)
      CALL ABRT
      STOP
 9008 FORMAT(31H ROUNDOFF ERROR IN BEND - STOP )
      END
C*MODULE ZMATRX  *DECK BMAT
      SUBROUTINE BMAT(B,IZMAT,ZMAT,S,T,NZVAR,NVAR,NCOORD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL OUT,GOPARR,DSKWRK,MASWRK
C
      DIMENSION B(NZVAR,NCOORD),IZMAT(*),ZMAT(*),S(NZVAR,NZVAR),T(NZVAR)
C
      PARAMETER (MXATM=2000)
C
      COMMON /INFOA / NAT,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 /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00)
      DATA DBUGME/8HBINV    /
C
      OUT = (NPRINT.EQ.1 .OR. EXETYP.EQ.DBUGME)  .AND. MASWRK
C
C     ----- CONSTRUCT THE B MATRIX -----
C     A WORD TO THE WISE, THE INITIAL B MATRIX IS NZVAR X NCOORD,
C     BUT IT ENDS UP BEING NVAR X NCOORD, IF NVAR.LT.NZVAR
C
      CALL VCLR(B,1,3*NAT*NZVAR)
C
C     ----- CONSTRUCT THE B MATRIX -----
C
      NLIBE=0
      NOINT = 0
      IADD = 1
C
  110 CONTINUE
      NOINT=NOINT+1
      IF (NOINT.GT.NZVAR) GO TO 230
C
      ITYPE = IZMAT(IADD)
      GO TO (120,140,160,180,200,210,215), ITYPE
C
C     ----- BOND STRETCH -----
C
  120 CONTINUE
      CALL BSTR(ZMAT(NOINT),NOINT,IZMAT(IADD+1),
     *          IZMAT(IADD+2),C,B,NZVAR)
      IADD = IADD + 3
C
      GO TO 110
C
C     ----- ANGLE BEND -----
C
  140 CONTINUE
      CALL BEND(ZMAT(NOINT),NOINT,IZMAT(IADD+1),IZMAT(IADD+2),
     *          IZMAT(IADD+3),C,B,NZVAR)
      IADD = IADD + 4
      GO TO 110
C
C     ----- TORSION -----
C     DIHEDRAL ANGLE BETWEEN TWO PLANES HAVING TWO ATOMS IN COMMON
C
  160 CONTINUE
      CALL TORS(ZMAT(NOINT),NOINT,IZMAT(IADD+1),IZMAT(IADD+2),
     *          IZMAT(IADD+3),IZMAT(IADD+4),C,B,NZVAR)
      IADD = IADD + 5
      GO TO 110
C
C     ----- OUT OF PLANE ANGLE BEND -----
C
  180 CONTINUE
      CALL OPLA(ZMAT(NOINT),NOINT,IZMAT(IADD+1),IZMAT(IADD+2),
     *          IZMAT(IADD+3),IZMAT(IADD+4),C,B,NZVAR)
      IADD = IADD + 5
      GO TO 110
C
C     ----- LINEAR ANGLE BEND ----- (THIS IS 2 COORDINATES)
C
  200 CONTINUE
      NLIBE=NLIBE+1
      CALL LIBE(ZMAT(NOINT),NOINT,NLIBE,IZMAT(IADD+1),
     *          IZMAT(IADD+2),IZMAT(IADD+3),C,B,NZVAR)
      ZMAT(NOINT+1)=ZMAT(NOINT)
      NOINT=NOINT+1
      IADD = IADD+4
      GO TO 110
C
C     ----- DIHEDRAL ANGLE BETWEEN TWO PLANES HAVING 1 ATOM IN COMMON
C
  210 CONTINUE
      CALL BDIHA5(IW,NOINT,NZVAR,IZMAT(IADD+1),C,ZMAT(NOINT),B)
      IADD=IADD+6
      GO TO 110
C
C     ----- "GHOST ATOM" TORSION -----  (I1,I2,J,K,L1,L2)
C     DIHEDRAL ANGLE BETWEEN TWO PLANES HAVING A BOND IN COMMON
C     (FIRST PLANE DEFINED BY MEDIAN OF I1-J-I2 ANGLE AND J-K
C     BOND, SECOND PLANE DEFINED BY MEDIAN OF L1-K-L2 ANGLE
C     AND J-K BOND.
C
  215 CONTINUE
      CALL GTORS(ZMAT(NOINT),NOINT,IZMAT(IADD+1),IZMAT(IADD+2),
     *           IZMAT(IADD+3),IZMAT(IADD+4),IZMAT(IADD+5),
     *           IZMAT(IADD+6),C,B,NZVAR)
      IADD = IADD + 7
      GO TO 110
C
C               UNITS ARE BOHRS AND RADIANS
C
  230 CONTINUE
      CALL DAWRIT(IDAF,IODA,ZMAT,NZVAR      ,39,0)
      CALL DAWRIT(IDAF,IODA,B   ,3*NAT*NZVAR,41,0)
C
      IF(OUT) THEN
         DUM=ZERO
         CALL PRZDAT(NZVAR,IZMAT,ZMAT,DUM,0)
         WRITE (IW,9028)
         CALL PRSQ(B,3*NAT,NZVAR,NZVAR)
      END IF
C
C     ----- SYMMETRIZATION OF THE B MATRIX -----
C
      DO 260 I = 1,3*NAT
         DO 250 J = 1,NVAR
            DUM = ZERO
            DO 240 K = 1,NZVAR
               DUM = DUM + B(K,I)*S(K,J)
  240       CONTINUE
            T(J) = DUM
  250    CONTINUE
         DO 255 J = 1,NVAR
            B(J,I) =T(J)
  255    CONTINUE
  260 CONTINUE
C
C     ------ IF THE MATRIX CHANGES DIMENSIONALITY (NZVAR > NVAR)
C            WE "SQUEEZE" THE BMATRIX DOWN TO A SMALLER DIMENSION -----
C
      CALL SQZMAT(B,3*NAT,NVAR,NZVAR)
C
      IF (OUT) THEN
         WRITE(IW,9048)
         CALL PRSQ(B,3*NAT,NVAR,NVAR)
      END IF
C
C     ----- SYMMETRIZE THE INTERNAL COORDINATES -----
C
      DO 350 J = 1,NVAR
         T(J) = ZERO
         DO 340 K = 1,NZVAR
            T(J) = T(J) + ZMAT(K)*S(K,J)
  340    CONTINUE
  350 CONTINUE
      DO 355 J = 1,NVAR
         ZMAT(J) = T(J)
  355 CONTINUE
C
C
C     ----- SAVE THE SYMMETRIZED-B AND Z MATRICES ON DIRECT ACCESS -----
C
C
      CALL DAWRIT(IDAF,IODA,ZMAT,NVAR       ,47,0)
      CALL DAWRIT(IDAF,IODA,B   ,NVAR*NCOORD,49,0)
      RETURN
C
 9028 FORMAT(/10H B MATRIX /)
 9048 FORMAT(/,' SYMMETRIZED -B- MATRIX -----')
      END
C*MODULE ZMATRX  *DECK BSTR
      SUBROUTINE BSTR(EQVAL,NOINT,I,J,C,B,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     -----THIS ROUTINE COMPUTES THE B MATRIX ELEMENTS FOR A
C          BOND STRETCH AS DEFINED BY WILSON (SEE WDC P.55) -----
C
      DIMENSION C(3,1),B(NDIM,1)
      DIMENSION RIJ(3)
      DIMENSION CI(6,6),X(6),Y(6)
C
      PARAMETER (ZERO=0.0D+00)
C
      DIJSQ = ZERO
      DO 100 M = 1,3
         RIJ(M) = C(M,J)-C(M,I)
         DIJSQ = DIJSQ+RIJ(M)*RIJ(M)
  100 CONTINUE
      R  = SQRT(DIJSQ)
      DO 120 M = 1,3
         NOCOL1 = 3*(I-1)+M
         NOCOL2 = 3*(J-1)+M
         B(NOINT,NOCOL1) = -RIJ(M)/R
         B(NOINT,NOCOL2) = -B(NOINT,NOCOL1)
  120 CONTINUE
C
C  ADD ON QUADRATIC TERMS
C
      CALL CSTR(CI,C,J,I)
C
      X(1) = C(1,J)
      X(2) = C(2,J)
      X(3) = C(3,J)
      X(4) = C(1,I)
      X(5) = C(2,I)
      X(6) = C(3,I)
C
      CALL TFSQU(ANS,CI,X,Y,6,1)
      EQVAL = R + 0.5D+00*ANS
C
      RETURN
      END
C*MODULE ZMATRX  *DECK CSTR
      SUBROUTINE CSTR(CI,C,J,I)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION CI(6,6),C(3,*),RIJ(3)
C
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
C
C   QUADRATIC TERMS FOR STRETCHES
C
      DIJSQ = ZERO
      DO 100 M = 1,3
         RIJ(M) = C(M,J)-C(M,I)
         DIJSQ = DIJSQ+RIJ(M)*RIJ(M)
  100 CONTINUE
      R  = SQRT(DIJSQ)
      RIN = ONE/R
      CUBRIN = RIN**3
      DO 120 M = 1,3
         TEMP        =  RIN - (RIJ(M)*RIJ(M)*CUBRIN)
         CI(M,M)     =  TEMP
         CI(M+3,M+3) =  TEMP
         CI(M,M+3)   = -TEMP
         CI(M+3,M)   = -TEMP
  120 CONTINUE
C
C     CROSS TERMS
C     X-Y
C
      TEMP   = (-CUBRIN)*(C(1,J)-C(1,I))*(C(2,J)-C(2,I))
      CI(1,2) =  TEMP
      CI(4,5) =  TEMP
      CI(2,1) =  TEMP
      CI(5,4) =  TEMP
      CI(1,5) = -TEMP
      CI(4,2) = -TEMP
      CI(2,4) = -TEMP
      CI(5,1) = -TEMP
C
C   X-Z
C
      TEMP   = (-CUBRIN)*(C(1,J)-C(1,I))*(C(3,J)-C(3,I))
      CI(1,3) =  TEMP
      CI(4,6) =  TEMP
      CI(3,1) =  TEMP
      CI(6,4) =  TEMP
      CI(1,6) = -TEMP
      CI(4,3) = -TEMP
      CI(3,4) = -TEMP
      CI(6,1) = -TEMP
C
C   Y-Z
C
      TEMP   = (-CUBRIN)*(C(2,J)-C(2,I))*(C(3,J)-C(3,I))
      CI(2,3) =  TEMP
      CI(5,6) =  TEMP
      CI(3,2) =  TEMP
      CI(6,5) =  TEMP
      CI(2,6) = -TEMP
      CI(5,3) = -TEMP
      CI(3,5) = -TEMP
      CI(6,2) = -TEMP
C
      RETURN
      END
C*MODULE ZMATRX  *DECK FILLZM
      SUBROUTINE FILLZM(IX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      DIMENSION IX(*)
      PARAMETER (MXATM=2000)
      COMMON /CONNEC/ IZ(MXATM,4),IATMOR(MXATM),NATOMS,NRATM,NDATM
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- FILL Z-MATRIX USING $DATA INTERNAL DEFINITIONS -----
C
      IF(NDATM.GT.0) THEN
         IF (MASWRK) WRITE(IW,9010)
         CALL ABRT
         STOP
      END IF
C
      I = 1
      ICNT = 0
C
   20 CONTINUE
      ICNT = ICNT + 1
      IF(ICNT.GT.NATOMS) GO TO 100
      DO 30 N = 1,NATOMS
         IF(IATMOR(N).EQ.ICNT) GO TO 40
   30 CONTINUE
   40 CONTINUE
         IF(IATMOR(N).EQ.1) GO TO 20
         IX(I)    = 1
         IX(I+1)  = IATMOR(IATMOR(N))
         IX(I+2)  = IATMOR(IZ(IATMOR(N),1))
         IF(IATMOR(N).EQ.2) THEN
           I = I + 3
           GO TO 20
         END IF
         IX(I+3)  = 2
         IX(I+4)  = IATMOR(IATMOR(N))
         IX(I+5)  = IATMOR(IZ(IATMOR(N),1))
         IX(I+6)  = IATMOR(IZ(IATMOR(N),2))
         IF(IATMOR(N).EQ.3) THEN
           I = I + 7
           GO TO 20
         END IF
         ITEST = IABS(IZ(IATMOR(N),4))
         IF(ITEST.EQ.1) GO TO 50
         IF(ITEST.GE.2) CALL ABRT
         IX(I+7)  = 3
         IX(I+8)  = IATMOR(IATMOR(N))
         IX(I+9)  = IATMOR(IZ(IATMOR(N),1))
         IX(I+10) = IATMOR(IZ(IATMOR(N),2))
         IX(I+11) = IATMOR(IZ(IATMOR(N),3))
         I = I + 12
         GO TO 20
   50    CONTINUE
         IX(I+7)  = 2
         IX(I+8)  = IATMOR(IATMOR(N))
         IX(I+9)  = IATMOR(IZ(IATMOR(N),1))
         IX(I+10) = IATMOR(IZ(IATMOR(N),3))
         I = I + 11
         GO TO 20
  100 CONTINUE
      RETURN
C
 9010 FORMAT(//1X,'*** ERROR!'/1X,'GAMESS DOES NOT SUPPORT AUTOMATIC',
     *          ' $ZMAT GENERATION WITH DUMMY ATOMS.'/
     *       1X,'1. REWRITE YOUR $DATA WITHOUT ANY DUMMY ATOMS, OR'/
     *       1X,'2. RUN WITH NZVAR=0 TO SEE WHICH REAL ATOMS ARE',
     *          ' GENERATED,'/
     *       1X,'   AND THEN READ IN AN EXPLICIT $ZMAT FOR THE REAL',
     *          ' ATOMS.')
      END
C*MODULE ZMATRX  *DECK GINVR
      SUBROUTINE GINVR(B,BINV,EIGVEC,G,RECIPM,EIGVAL,SCR,IA,GLIN,
     *                 NVAR,NCOORD,NAT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION B(NVAR,NCOORD),EIGVEC(NVAR,NVAR),G(NVAR,NVAR),
     *          RECIPM(NCOORD),EIGVAL(NVAR),SCR(NVAR,8),IA(NVAR),
     *          GLIN(*),BINV(NCOORD,NVAR)
      DIMENSION DET(2)
C
      LOGICAL OUT,DECOMP,GOPARR,DSKWRK,MASWRK,FIRST
C
      PARAMETER (MXATM=2000)
C
      SAVE FIRST
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MASSES/ ZMASS(MXATM)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /VBDCMP/ DECOMP,NROW,NDEG
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, FOUR=4.0D+00, TOL=1.0D-05)
C
      DATA DBUGME/8HGINVR   /
      DATA FIRST/.TRUE./
      DATA VSCF  /8HVSCF    /
C
      OUT = MASWRK  .AND.  (EXETYP.EQ.DBUGME  .OR.  NPRINT.EQ.1)
C
C     ----- CALCULATE G = B * (M**-1) * (B)'
C
      I0=0
      DO 100 I=1,NAT
         IF(ZMASS(I).NE.ZERO) THEN
            TEMP = ONE/ZMASS(I)
         ELSE
            TEMP = ONE
         END IF
         RECIPM(I0+1) = TEMP
         RECIPM(I0+2) = TEMP
         RECIPM(I0+3) = TEMP
         I0=I0+3
  100 CONTINUE
C
      DO 130 I=1,NVAR
         DO 120 J=1,NVAR
            SUM = ZERO
            DO 110 K=1,NCOORD
            SUM = SUM + B(I,K)*RECIPM(K)*B(J,K)
  110       CONTINUE
            G(I,J) = SUM
  120    CONTINUE
  130 CONTINUE
C
      IF(OUT) THEN
         WRITE(IW,*) 'GINVR: G MATRIX IS'
         CALL PRSQ(G,NVAR,NVAR,NVAR)
      END IF
C
      IF (.NOT.DECOMP) GO TO 260
C
C     ----- PACK G INTO TRIANGULAR STORAGE -----
C
      DO 150 I=1,NVAR
         IROW = (I*I-I)/2
         DO 140 J=1,I
            GLIN(IROW+J) = G(I,J)
  140    CONTINUE
  150 CONTINUE
C
C     ----- DIAGONALIZE G -----
C
      CALL GLDIAG(NVAR,NVAR,NVAR,GLIN,SCR,EIGVAL,EIGVEC,IERR,IA)
C
C     ----- COUNT NO. OF ZERO EIGENVALUES, FIND G**-1, AND PUT IT IN G.
C
      NZERO = 0
      DO 160 I=1,NVAR
         IF (ABS(EIGVAL(I)) .LT. TOL) NZERO = NZERO + 1
  160 CONTINUE
C
      I3NM6 = NCOORD - NDEG
      INDEPC = NVAR - NZERO
      IF (INDEPC .LT. I3NM6) THEN
         IF (MASWRK) WRITE(IW,9077) INDEPC,I3NM6
         CALL ABRT
      END IF
C
      NZERO = NZERO + 1
      DO 220 I=1,NVAR
         DO 210 J=1,NVAR
            SUM = ZERO
            DO 200 K=NZERO,NVAR
               SUM = SUM + EIGVEC(I,K)*EIGVEC(J,K)/EIGVAL(K)
  200       CONTINUE
            G(I,J) = SUM
  210    CONTINUE
  220 CONTINUE
      GO TO 290
C
C     ----- CONSTRUCT THE G INVERSE MATRIX -----
C           WHEN NOT DOING A DECOMPOSITION
C
  260 CONTINUE
      INFO=0
      CALL DGEFA(G,NVAR,NVAR,IA,INFO)
      IF(INFO.NE.0) THEN
         IF (MASWRK) THEN
         WRITE(IW,*) 'YOUR B-MATRIX IS SINGULAR, BOMBING...'
         WRITE(IW,*) 'YOU SHOULD CHECK YOUR Z-MATRIX INPUT CAREFULLY.'
         END IF
         CALL ABRT
         STOP
      END IF
      CALL DGEDI(G,NVAR,NVAR,IA,DET,SCR,11)
C
C     ---- TEST THE MAGNITUDE OF THE DETERMINANT FOR LINEAR COORDINATE
C          QUALITY [DET(2) IS THE POWER OF TEN] ----
C
      IDET = INT( DET(2) )
      RANGE = DET(2)/NVAR
      IF (FIRST) THEN
         IF(MASWRK) WRITE(IW,9080) IDET
         FIRST=.FALSE.
      END IF
      IF (RANGE.LT.-ONE  .AND.  RANGE.GT.-FOUR) THEN
         IF(MASWRK) WRITE(IW,9080) IDET
         IF(MASWRK) WRITE(IW,9090)
      ELSE IF (RANGE.LT.-FOUR) THEN
         IF(MASWRK) WRITE(IW,9080) IDET
         IF(MASWRK) WRITE(IW,9100)
         IF(RUNTYP.EQ.VSCF) RETURN
         CALL ABRT
      END IF
C
C     ---- FINALLY, CALCULATE B**-1' = (G**-1) * B * (M**-1)
C
  290 CONTINUE
      IF(OUT) THEN
         WRITE(IW,*) 'GINVR: G INVERSE MATRIX IS'
         CALL PRSQ(G,NVAR,NVAR,NVAR)
      END IF
C
      DO 320 I = 1,NVAR
         DO 310 J = 1,NCOORD
            SUM = ZERO
            DO 300 K = 1,NVAR
               SUM = SUM + G(I,K)*B(K,J)
  300       CONTINUE
           BINV(J,I) = SUM*RECIPM(J)
  310    CONTINUE
  320 CONTINUE
C
      CALL DAWRIT(IDAF,IODA,BINV,NCOORD*NVAR,50,0)
      IF(OUT) THEN
         WRITE(IW,*) 'GINVR: B INVERSE MATRIX IS'
         CALL PRSQ(BINV,NVAR,NCOORD,NCOORD)
      END IF
      RETURN
C
 9077 FORMAT(' **** ERROR IN GINVR *** '/
     *       ' THERE ARE ONLY ',I3,' INDEPENDENT INTERNAL COORDINATES.'/
     *       ' YOU MUST PROVIDE ',I3)
 9080 FORMAT(/1X,'THE DETERMINANT OF THE G MATRIX IS 10**(',I6,')')
 9090 FORMAT(1X,'* * * WARNING * * *'/
     *       1X,'YOUR B MATRIX HAS A SMALL DETERMINANT.'/
     *       1X,'THIS SUGGESTS YOU MAY HAVE BADLY CHOSEN COORDINATES.')
 9100 FORMAT(1X,'SIZE OF DETERMINANT SUGGESTS A SINGULAR G MATRIX.'/
     *       1X,'BOMBING.  CHECK YOUR COORDINATES')
      END
C*MODULE ZMATRX  *DECK GTBOND
      SUBROUTINE GTBOND(IPAIRS,MAXBND,NBONDS,LUNIT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION MAKEFP
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      CHARACTER*8 FRGNME
C
      DIMENSION IPAIRS(2,MAXBND)
      DIMENSION RCOV(86)
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXPT=2000, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FRGMSS/ FPMASS(MXPT),FMASS(MXFGPT),
     *                FPNUC(MXPT),FGNUC(MXFGPT)
      COMMON /INFOA / NAT,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 /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      PARAMETER (ZERO=0.0D+00, TOANGS=0.52917724924D+00)
C
C      COVALENT RADII FROM J.EMSLEY, "THE ELEMENTS", 2ND EDITION, 1991
C      EXCEPT VAN DER WAALS RADII FOR HE,NE,AR,KR (SAME SOURCE),
C      AND GUESSES FOR NA,V,CR,RB,TC,PM,EU,YB,AT,RN
C
      DATA (RCOV(NUCZ),NUCZ=1,2)/0.30D+00,1.22D+00/
      DATA (RCOV(NUCZ),NUCZ=3,10)
     *  /1.23D+00,0.89D+00,0.88D+00,0.77D+00,
     *   0.70D+00,0.66D+00,0.58D+00,1.60D+00/
      DATA (RCOV(NUCZ),NUCZ=11,18)
     *  /1.66D+00,1.36D+00,1.25D+00,1.17D+00,
     *   1.10D+00,1.04D+00,0.99D+00,1.91D+00/
      DATA (RCOV(NUCZ),NUCZ=19,36)
     *  /2.03D+00,1.74D+00,
     *   1.44D+00,1.32D+00,1.22D+00,1.19D+00,1.17D+00,
     *   1.165D+00,1.16D+00,1.15D+00,1.17D+00,1.25D+00,
     *   1.25D+00,1.22D+00,1.21D+00,1.17D+00,1.14D+00,1.98D+00/
      DATA (RCOV(NUCZ),NUCZ=37,54)
     *  /2.22D+00,1.92D+00,
     *   1.62D+00,1.45D+00,1.34D+00,1.29D+00,1.27D+00,
     *   1.24D+00,1.25D+00,1.28D+00,1.34D+00,1.41D+00,
     *   1.50D+00,1.40D+00,1.41D+00,1.37D+00,1.33D+00,2.09D+00/
      DATA (RCOV(NUCZ),NUCZ=55,86)
     *  /2.35D+00,1.98D+00,
     *   1.69D+00,1.65D+00,1.65D+00,1.64D+00,1.65D+00,1.66D+00,1.65D+00,
     *   1.61D+00,1.59D+00,1.59D+00,1.58D+00,1.57D+00,1.56D+00,1.56D+00,
     *   1.56D+00,1.44D+00,1.34D+00,1.30D+00,1.28D+00,
     *   1.26D+00,1.26D+00,1.29D+00,1.34D+00,1.44D+00,
     *   1.55D+00,1.54D+00,1.52D+00,1.53D+00,1.50D+00,2.20D+00/
C
      DATA MAKEFP/8HMAKEFP  /
C
C     ----- THREE OPTIONS, CONTROLLED BY -LUNIT- -----
C     LUNIT=-1: JUST COUNT THE NUMBER OF BONDS AND RETURN IN -NBONDS-,
C               WITHOUT STORING THEM IN THE BONDED ATOMS ARRAYS.
C               STORAGE FOR -IPAIRS- NEED NOT BE ALLOCATED.
C     LUNIT=0:  FILL -IPAIRS- WITH LIST OF BONDED ATOMS AND
C               RETURN -NBONDS-
C     LUNIT>0:  WRITE THE PREVIOUSLY DEFINED LIST OF BONDED ATOMS
C               TO -LUNIT-
C
      IF(LUNIT.GT.0) GO TO 600
C
C     ----- FILL -IPAIRS- -----
C
      NBONDS = 0
      DO 130 I=1,NAT
         INUCZ = INT(ZAN(I)) + IZCORE(I)
         IF(RUNTYP.EQ.MAKEFP  .AND.  INUCZ.EQ.0) GO TO 130
         RADI = 1.6D+00
         IF(INUCZ.EQ.1)                 RADI =         RCOV(1)
         IF(INUCZ.GT.1.AND.INUCZ.LE.86) RADI = 1.2D+00*RCOV(INUCZ)
         DO 120 J=1,I-1
            JNUCZ = INT(ZAN(J)) + IZCORE(J)
            IF(RUNTYP.EQ.MAKEFP  .AND.  JNUCZ.EQ.0) GO TO 130
            RADJ = 1.6D+00
            IF(JNUCZ.EQ.1)                 RADJ =         RCOV(1)
            IF(JNUCZ.GT.1.AND.JNUCZ.LE.86) RADJ = 1.2D+00*RCOV(JNUCZ)
            DIST = SQRT((C(1,I)-C(1,J))**2
     *                + (C(2,I)-C(2,J))**2
     *                + (C(3,I)-C(3,J))**2) * TOANGS
            BOND = RADI + RADJ
            IF(DIST.LE.BOND) THEN
               NBONDS = NBONDS + 1
               IF(LUNIT.EQ.0) THEN
                  IF(NBONDS.GT.MAXBND) THEN
                     IF(MASWRK) WRITE(IW,9000) MAXBND
                     CALL ABRT
                  END IF
                  IPAIRS(1,NBONDS) = I
                  IPAIRS(2,NBONDS) = J
               END IF
            END IF
  120    CONTINUE
  130 CONTINUE
C
      IF(NFRG.EQ.0) RETURN
C                   ******
C
      IAT  = 0
      IIAT = 0
      JJAT = 0
      JAT0 = 0
      JJAT0= 0
      KMASS=0
      DO 240 IFRG=1,NFRG
         DO 220 III=1,NMPTS(IFRG)
            IAT=IAT+1
            IF(FMASS(III+KMASS).LE.ZERO) GO TO 220
            IIAT = IIAT+1
            INUCZ = INT(FGNUC(III+KMASS)+0.001D+00)
            RADI = 1.6D+00
            IF(INUCZ.EQ.1) RADI = RCOV(1)
            IF(INUCZ.GT.1.AND.INUCZ.LE.86) RADI = 1.2D+00*RCOV(INUCZ)
            JAT = JAT0
            JJAT= JJAT0
            DO 210 JJJ=1,III-1
               JAT=JAT+1
               IF(FMASS(JJJ+KMASS).LE.ZERO) GO TO 210
               JJAT = JJAT+1
               JNUCZ = INT(FGNUC(JJJ+KMASS)+0.001D+00)
               RADJ = 1.6D+00
               IF(JNUCZ.EQ.1) RADJ = RCOV(1)
               IF(JNUCZ.GT.1.AND.JNUCZ.LE.86) RADJ = 1.2D+00*RCOV(JNUCZ)
               DIST = SQRT((EFC(1,IAT)-EFC(1,JAT))**2
     *                   + (EFC(2,IAT)-EFC(2,JAT))**2
     *                   + (EFC(3,IAT)-EFC(3,JAT))**2) * TOANGS
               BOND = RADI + RADJ
               IF(DIST.LE.BOND) THEN
                  NBONDS = NBONDS + 1
                  IF(LUNIT.EQ.0) THEN
                     IF(NBONDS.GT.MAXBND) THEN
                        IF(MASWRK) WRITE(IW,9000) MAXBND
                        CALL ABRT
                     END IF
                     IPAIRS(1,NBONDS) = IIAT+NAT
                     IPAIRS(2,NBONDS) = JJAT+NAT
                  END IF
               END IF
  210       CONTINUE
  220    CONTINUE
         JAT0 = IAT
         JJAT0=IIAT
         KMASS=KMASS + NMPTS(IFRG)
  240 CONTINUE
      RETURN
C     ******
C
C     ----- PUNCH LIST OF BONDED ATOMS -----
C
  600 CONTINUE
      IF(.NOT.MASWRK) RETURN
      IF(NBONDS.LE.6) THEN
         WRITE(LUNIT,8010) (IPAIRS(1,I),IPAIRS(2,I),I=1,NBONDS)
      ELSE
         WRITE(LUNIT,8020) (IPAIRS(1,I),IPAIRS(2,I),I=1,6)
         WRITE(LUNIT,8030) (IPAIRS(1,I),IPAIRS(2,I),I=7,NBONDS)
      END IF
      RETURN
C
 8010 FORMAT('BONDATOMS ',6(I4,I4,2X))
 8020 FORMAT('BONDATOMS ',6(I4,I4,2X),' >')
 8030 FORMAT(7(I4,I4,2X),:,' >')
 9000 FORMAT(1X,'GTBOND: NOT ENOUGH BOND STORAGE ALLOCATED:',I10)
      END
C*MODULE ZMATRX  *DECK GTORS
      SUBROUTINE GTORS(EQVAL,NOINT,I1,I2,J,K,L1,L2,C,B,NDIM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     -----THIS ROUTINE COMPUTES THE B MATRIX ELEMENTS FOR A
C          "GHOST ATOM" TORSION, SIMILAR TO THE TORSION DEF. BY WILSON.
C          SEE PP. 60-61 OF "MOLECULAR VIBRATIONS" BY E.B. WILSON,
C          J.C. DECIUS AND P.C. CROSS, MCGRAW-HILL, NY, 1955.
C
C     -----I1 AND I2 ARE THE END ATOMS OF THE I1-J-I2 PLANE.  I2 IS
C          NEARER THE OBSERVER.  L1, L2 ARE THE END ATOMS OF THE L1-K-L2
C          PLANE, L2 IS NEARER THE OBSERVER.  THE TORSION ANGLE IS THE
C          ANGLE BETWEEN THE PLANE DEFINED BY THE MEDIAN OF THE I1-J-I2
C          ANGLE AND THE J-K BOND, AND THE PLANE DEFINED BY THE MEDIAN
C          OF THE L1-K-L2 ANGLE AND THE J-K BOND.
C
      LOGICAL I1EQI2,L1EQL2,GOPARR,DSKWRK,MASWRK
C
      DIMENSION C(3,1),B(NDIM,1)
      DIMENSION RIJ(3),RJK(3),RKL(3),EIJ(3),EJK(3),EKL(3),CR1(3),
     + CR2(3)
      DIMENSION RI1J(3),RI2J(3),EI1J(3),EI2J(3),CR3(3),CR4(3),
     +          RKL1(3),RKL2(3),EKL1(3),EKL2(3),CR5(3),CR6(3)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, PT5=0.5D+00)
      PARAMETER (TOLRD=1.0001D+00, TOL=1.0D-07)
C
C
      PI = ACOS(-ONE)
      DIJSQ  = ZERO
      DJKSQ  = ZERO
      DKLSQ  = ZERO
      DI1JSQ = ZERO
      DI2JSQ = ZERO
      DKL1SQ = ZERO
      DKL2SQ = ZERO
      I1EQI2 = I1 .EQ. I2
      L1EQL2 = L1 .EQ. L2
      DO 120 M = 1,3
         RIJ(M)  = C(M,J) - PT5*(C(M,I1) + C(M,I2))
         DIJSQ   =    DIJSQ + RIJ(M)*RIJ(M)
         RJK(M)  = C(M,K) - C(M,J)
         DJKSQ   =    DJKSQ + RJK(M)*RJK(M)
         RKL(M)  = PT5*(C(M,L1) + C(M,L2)) - C(M,K)
         DKLSQ   =    DKLSQ + RKL(M)*RKL(M)
         RI1J(M) = C(M,J) - C(M,I1)
         DI1JSQ  =    DI1JSQ + RI1J(M)*RI1J(M)
         RI2J(M) = C(M,J) - C(M,I2)
         DI2JSQ  =    DI2JSQ + RI2J(M)*RI2J(M)
         RKL1(M) = C(M,L1) - C(M,K)
         DKL1SQ  =    DKL1SQ + RKL1(M)*RKL1(M)
         RKL2(M) = C(M,L2) - C(M,K)
         DKL2SQ  =    DKL2SQ + RKL2(M)*RKL2(M)
  120 CONTINUE
      DIJ  = SQRT(DIJSQ)
      DJK  = SQRT(DJKSQ)
      DKL  = SQRT(DKLSQ)
      DI1J = SQRT(DI1JSQ)
      DI2J = SQRT(DI2JSQ)
      DKL1 = SQRT(DKL1SQ)
      DKL2 = SQRT(DKL2SQ)
      DO 180 M = 1,3
         EIJ(M)  = RIJ(M)/DIJ
         EJK(M)  = RJK(M)/DJK
         EKL(M)  = RKL(M)/DKL
         EI1J(M) = RI1J(M)/DI1J
         EI2J(M) = RI2J(M)/DI2J
         EKL1(M) = RKL1(M)/DKL1
         EKL2(M) = RKL2(M)/DKL2
  180 CONTINUE
C
      CR1(1) = EIJ(2)*EJK(3)-EIJ(3)*EJK(2)
      CR1(2) = EIJ(3)*EJK(1)-EIJ(1)*EJK(3)
      CR1(3) = EIJ(1)*EJK(2)-EIJ(2)*EJK(1)
C
      CR2(1) = EJK(2)*EKL(3)-EJK(3)*EKL(2)
      CR2(2) = EJK(3)*EKL(1)-EJK(1)*EKL(3)
      CR2(3) = EJK(1)*EKL(2)-EJK(2)*EKL(1)
C
      CR3(1) = EI1J(2)*EJK(3)-EI1J(3)*EJK(2)
      CR3(2) = EI1J(3)*EJK(1)-EI1J(1)*EJK(3)
      CR3(3) = EI1J(1)*EJK(2)-EI1J(2)*EJK(1)
C
      CR4(1) = EI2J(2)*EJK(3)-EI2J(3)*EJK(2)
      CR4(2) = EI2J(3)*EJK(1)-EI2J(1)*EJK(3)
      CR4(3) = EI2J(1)*EJK(2)-EI2J(2)*EJK(1)
C
      CR5(1) = EJK(2)*EKL1(3)-EJK(3)*EKL1(2)
      CR5(2) = EJK(3)*EKL1(1)-EJK(1)*EKL1(3)
      CR5(3) = EJK(1)*EKL1(2)-EJK(2)*EKL1(1)
C
      CR6(1) = EJK(2)*EKL2(3)-EJK(3)*EKL2(2)
      CR6(2) = EJK(3)*EKL2(1)-EJK(1)*EKL2(3)
      CR6(3) = EJK(1)*EKL2(2)-EJK(2)*EKL2(1)
C
      DOTPJ  = ZERO
      DOTPK  = ZERO
      DOTI1J = ZERO
      DOTI2J = ZERO
      DOTKL1 = ZERO
      DOTKL2 = ZERO
      DO 220 M = 1,3
         DOTPJ  = DOTPJ-EIJ(M)*EJK(M)
         DOTPK  = DOTPK-EJK(M)*EKL(M)
         DOTI1J = DOTI1J-EI1J(M)*EJK(M)
         DOTI2J = DOTI2J-EI2J(M)*EJK(M)
         DOTKL1 = DOTKL1-EKL1(M)*EJK(M)
         DOTKL2 = DOTKL2-EKL2(M)*EJK(M)
  220 CONTINUE
      IF (ONE.LE.ABS(DOTPJ))  GO TO 320
      IF (ONE.LE.ABS(DOTPK))  GO TO 320
C
      IF (ONE.LT.ABS(DOTI1J)) GO TO 320
      IF (ONE.LT.ABS(DOTI2J)) GO TO 320
      IF (ONE.LT.ABS(DOTKL1)) GO TO 320
      IF (ONE.LT.ABS(DOTKL2)) GO TO 320
C
      SINPJ  = SQRT(ONE-DOTPJ*DOTPJ)
      SINPK  = SQRT(ONE-DOTPK*DOTPK)
      SINI1J = SQRT(ONE-DOTI1J*DOTI1J)
      SINI2J = SQRT(ONE-DOTI2J*DOTI2J)
      SINKL1 = SQRT(ONE-DOTKL1*DOTKL1)
      SINKL2 = SQRT(ONE-DOTKL2*DOTKL2)
C
      DO 280 M = 1,3
         SMI1 = -PT5*CR3(M)/(DI1J*SINI1J*SINI1J)
         NOCOL1 = 3*(I1-1)+M
         B(NOINT,NOCOL1) = SMI1
C
         SMI2 = -PT5*CR4(M)/(DI2J*SINI2J*SINI2J)
         IF (I1EQI2) SMI2 = SMI2+SMI2
         NOCOL1 = 3*(I2-1)+M
         B(NOINT,NOCOL1) = SMI2
C
         F1 = - CR3(M)*DOTI1J/(DJK*SINI1J*SINI1J)
     *        - CR4(M)*DOTI2J/(DJK*SINI2J*SINI2J)
         F2 = - CR5(M)*DOTKL1/(DJK*SINKL1*SINKL1)
     *        - CR6(M)*DOTKL2/(DJK*SINKL2*SINKL2)
         F1 = PT5*F1
         F2 = PT5*F2
         SMJ = F1 + F2 - (SMI1 + SMI2)
         IF (I1EQI2) SMJ = F1 + F2 - SMI2
         NOCOL2 = 3*(J-1)+M
         B(NOINT,NOCOL2) = SMJ
C
         SML1 = PT5*CR5(M)/(DKL1*SINKL1*SINKL1)
         NOCOL3 = 3*(L1-1)+M
         B(NOINT,NOCOL3) = SML1
C
         SML2 = PT5*CR6(M)/(DKL2*SINKL2*SINKL2)
         NOCOL3 = 3*(L2-1)+M
         IF (L1EQL2) SML2 = SML2+SML2
         B(NOINT,NOCOL3) = SML2
C
         NOCOL4 = 3*(K-1)+M
         SUMNEG = -(SMI1+SMI2+SMJ+SML1+SML2)
         IF (I1EQI2) SUMNEG = SUMNEG+SMI1
         IF (L1EQL2) SUMNEG = SUMNEG+SML1
         B(NOINT,NOCOL4) = SUMNEG
  280 CONTINUE
      DOT = ZERO
      DO 300 M = 1,3
  300 DOT = DOT+(CR1(M)*CR2(M))/(SINPJ*SINPK)
      IF (ABS(DOT) .GT. TOLRD) GO TO 320
         IF (ABS(DOT) .GT. ONE) DOT = SIGN(ONE,DOT)
         EQVAL =   ACOS(DOT)
         IF(ABS(EQVAL).LT.TOL) EQVAL=ZERO
         IF(ABS(EQVAL-PI).LT.TOL) EQVAL=PI
         RETURN
C
  320 CONTINUE
      IF (MASWRK) WRITE (IW,9008)
      CALL ABRT
      STOP
 9008 FORMAT(31H ROUNDOFF ERROR IN GTORS - STOP )
      END
C*MODULE ZMATRX  *DECK LIBE
      SUBROUTINE LIBE(EQVAL,NOINT,NLIBE,I,J,K,C,B,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      PARAMETER (TOL1=1.0D-04, TOL2=5.0D-04)
      PARAMETER (MAXBE=15, NNAM=1)
C
      DIMENSION B(NDIM,1),C(3,2)
      DIMENSION QNAM(NNAM),KQNAM(NNAM),SLIBE(1)
      DIMENSION A(3),RJI(3),RJK(3),UN(3),APTS(3*MAXBE),RJA(3),
     *          UNIT(3),UP(3),EJI(3),EJK(3)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      SAVE APTS
C
      DATA SLIBE/8HLIBE    /, QNAM/8HAPTS    /
C
C     ----- THIS ROUTINE COMPUTES THE B MATRIX ELEMENTS FOR ONE OR
C           BOTH OF 2 PERPENDICULAR LINEAR BENDING COORDINATES.
C           SEE FERIGLE AND MEISTER, J CHEM PHYS 19,982(1951) -----
C
C     ----- I AND K ARE THE END ATOMS AND J IS THE CENTRAL ATOM. A
C           GIVES THE CARTESIAN COORDINATES OF A POINT IN SPACE SUCH
C           THAT THE VECTOR FROM ATOM J TO POINT A IS PERPENDICULAR TO
C           THE LINE I-J-K AND SERVES TO ORIENT THE COORDS IN SPACE.
C           THE FIRST OF THE TWO INTERNAL COORDINATES IS IN THE I-A-K
C           PLANE AND THE SECOND IN A PLANE PERPENDICULAR TO THE FIRST,
C           THROUGH POINTS I,J, AND K -----
C
      PI = ACOS(-ONE)
C
C     ----- READ THE PLANE CONTAINING THE POINTS A -----
C
      IF(NLIBE.GT.1) GO TO 100
         KQNAM(1) = 3*10*MAXBE + 3
         JRET=0
         CALL NAMEIO(IR,JRET,SLIBE,NNAM,QNAM,KQNAM,
     *               APTS,
     *               0,0,0,
     *    0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *    0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *    0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
         IF (JRET .NE. 1) GO TO 100
            IF (MASWRK) WRITE(IW,9000)
            CALL ABRT
  100 CONTINUE
C
      IF(NLIBE.GT.MAXBE) THEN
         IF (MASWRK) WRITE(IW,9030) MAXBE
         CALL ABRT
      END IF
C
C     ----- TAKE THE POINT A LYING IN THE PLANE, AND GENERATE
C           THE POINT A SUCH THAT A-J IS PERP. TO I-J-K -----
C
      DJKSQ=ZERO
      DOT=ZERO
      LOC=3*NLIBE-3
      DO 120 M=1,3
         A(M) = APTS(LOC+M)
         RJK(M) = C(M,K) - C(M,J)
         DJKSQ = DJKSQ + RJK(M)*RJK(M)
         RJA(M) = A(M) - C(M,J)
         DOT = DOT + RJK(M)*RJA(M)
  120 CONTINUE
      FACT=DOT/DJKSQ
      RJASQ = ZERO
      DO 140 M=1,3
         RJA(M) = RJA(M) - FACT*RJK(M)
         RJASQ = RJASQ + RJA(M)*RJA(M)
         A(M) = C(M,J) + RJA(M)
  140 CONTINUE
      IF(ABS(RJASQ).LT.TOL1) THEN
         IF (MASWRK) WRITE(IW,9010) (APTS(LOC+M),M=1,3)
         CALL ABRT
      END IF
C
C     ----- A IS NOW CORRECTLY LOCATED -----
C
      DJISQ = ZERO
      DJKSQ = ZERO
      DAJSQ = ZERO
      DO 210 M = 1,3
         RJI(M) = C(M,I)-C(M,J)
         DJISQ = DJISQ+RJI(M)*RJI(M)
         RJK(M) = C(M,K)-C(M,J)
         DJKSQ = DJKSQ+RJK(M)*RJK(M)
         UN(M) = A(M)-C(M,J)
         DAJSQ = DAJSQ+UN(M)*UN(M)
  210 CONTINUE
      DJI = SQRT(DJISQ)
      DJK = SQRT(DJKSQ)
      DAJ = SQRT(DAJSQ)
      DOTJ = ZERO
      DOTP = ZERO
      DO 230 M = 1,3
         EJI(M) = RJI(M)/DJI
         EJK(M) = RJK(M)/DJK
         UNIT(M) = UN(M)/DAJ
         DOTJ = DOTJ+EJI(M)*EJK(M)
         DOTP = DOTP+EJI(M)*UNIT(M)
  230 CONTINUE
      TEST = (ABS(DOTJ)-ONE)
      IF ((ABS(TEST).GT.TOL1) .OR.
     *    (ABS(DOTP).GT.TOL2)) THEN
         IF (MASWRK) WRITE (IW,9020)  TEST,DOTP
         CALL ABRT
      END IF
      IF (ABS(DOTJ)-ONE .GT. ZERO) DOTJ=SIGN(ONE,DOTJ)
      EQVAL = ACOS(DOTJ)
      IF(ABS(EQVAL)   .LT.TOL1) EQVAL=ZERO
      IF(ABS(EQVAL-PI).LT.TOL1) EQVAL=PI
C
C     ----- FIRST COMPONENT OF THE BEND -----
C
      NOCOL1 = 3*(I-1)
      NOCOL2 = 3*(K-1)
      NOCOL3 = 3*(J-1)
      DO 300 M = 1,3
         B(NOINT,NOCOL1+M) = -UNIT(M)/DJI
         B(NOINT,NOCOL2+M) = -UNIT(M)/DJK
         B(NOINT,NOCOL3+M) = (ONE/DJI+ONE/DJK)*UNIT(M)
  300 CONTINUE
C
C     ----- SECOND COMPONENT OF THE BEND -----
C
      UP(1) = EJK(2)*UNIT(3)-EJK(3)*UNIT(2)
      UP(2) = EJK(3)*UNIT(1)-EJK(1)*UNIT(3)
      UP(3) = EJK(1)*UNIT(2)-EJK(2)*UNIT(1)
      DO 400 M = 1,3
         B1 = -UP(M)/DJI
         B2 = -UP(M)/DJK
         B(NOINT+1,NOCOL1+M) = B1
         B(NOINT+1,NOCOL2+M) = B2
         B(NOINT+1,NOCOL3+M) = -(B1+B2)
  400 CONTINUE
      RETURN
C
 9000 FORMAT(1X,'NO $LIBE INPUT DATA FOUND...THIS IS REQUIRED')
 9010 FORMAT(1X,'THE POINT A IS COLINEAR WITH I,J,K'/
     *    1X,'A=',3F15.8)
 9020 FORMAT(1X,'*** ROUNDOFF ERROR IN LIBE - TEST,DOTP=',2E20.10)
 9030 FORMAT(1X,'ONLY',I4,' LINEAR BENDS MAY BE USED --- ERROR')
      END
C*MODULE ZMATRX  *DECK OPLA
      SUBROUTINE OPLA(EQVAL,NOINT,I,J,K,L,C,B,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
C     -----THIS ROUTINE COMPUTES THE B MATRIX ELEMENTS FOR AN
C          OUT OF PLANE WAGGING COORDINATE AS DEFINED BY WILSON.
C          SEE PP. 58-60 OF "MOLECULAR VIBRATIONS" BY E.B. WILSON,
C          J.C. DECIUS AND P.C. CROSS, MCGRAW-HILL, NY, 1955. ----
C
C     -----I IS THE END ATOM,, J IS THE APEX ATOM,, AND K AND L ARE
C         THE ANCHOR ATOMS -----
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DIMENSION  B(NDIM,1),C(3,1)
      DIMENSION RJI(3),RJK(3),RJL(3),EJI(3),EJK(3),EJL(3),
     +          C1(3),C2(3),C3(3)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TOL=1.0D-07)
C
      DATA CHECK/8HCHECK   /
C
      PIOVR2 = ZERO
      PIOVR2 = ACOS(PIOVR2)
      PI = PIOVR2+PIOVR2
      DJISQ = ZERO
      DJKSQ = ZERO
      DJLSQ = ZERO
      DO 100 M = 1,3
         RJI(M) = C(M,I)-C(M,J)
         DJISQ = DJISQ+RJI(M)*RJI(M)
         RJK(M) = C(M,K)-C(M,J)
         DJKSQ = DJKSQ+RJK(M)*RJK(M)
         RJL(M) = C(M,L)-C(M,J)
         DJLSQ = DJLSQ+RJL(M)*RJL(M)
  100 CONTINUE
      DJI = SQRT(DJISQ)
      DJK = SQRT(DJKSQ)
      DJL = SQRT(DJLSQ)
      DO 120 M = 1,3
         EJI(M) = RJI(M)/DJI
         EJK(M) = RJK(M)/DJK
         EJL(M) = RJL(M)/DJL
  120 CONTINUE
C
      C1(1) = EJK(2)*EJL(3)-EJK(3)*EJL(2)
      C1(2) = EJK(3)*EJL(1)-EJK(1)*EJL(3)
      C1(3) = EJK(1)*EJL(2)-EJK(2)*EJL(1)
C
      C2(1) = EJL(2)*EJI(3)-EJL(3)*EJI(2)
      C2(2) = EJL(3)*EJI(1)-EJL(1)*EJI(3)
      C2(3) = EJL(1)*EJI(2)-EJL(2)*EJI(1)
C
      C3(1) = EJI(2)*EJK(3)-EJI(3)*EJK(2)
      C3(2) = EJI(3)*EJK(1)-EJI(1)*EJK(3)
      C3(3) = EJI(1)*EJK(2)-EJI(2)*EJK(1)
C
      DOTI = ZERO
      DO 140 M = 1,3
  140 DOTI = DOTI+EJK(M)*EJL(M)
      IF (ONE.LE.ABS(DOTI)) GO TO 800
      SINI = SQRT(ONE-DOTI*DOTI)
      DET = EJI(1)*C1(1)+EJI(2)*C1(2)+EJI(3)*C1(3)
      SINT = DET/SINI
      IF (ONE.LE.ABS(SINT)) GO TO 810
      COST = SQRT(ONE-SINT*SINT)
      TANT = SINT/COST
      DO 200 M = 1,3
         NOCOL1 = 3*(I-1)+M
         SMI = ((C1(M)/(COST*SINI))-(TANT*EJI(M)))/DJI
         B(NOINT,NOCOL1) = SMI
         NOCOL2 = 3*(K-1)+M
         SMK = ((C2(M)/(COST*SINI))-((TANT*(EJK(M)-DOTI*EJL(M)))/(SINI*
     +        SINI)))/DJK
         B(NOINT,NOCOL2) = SMK
         NOCOL3 = 3*(L-1)+M
         SML = ((C3(M)/(COST*SINI))-((TANT*(EJL(M)-DOTI*EJK(M)))/(SINI*
     +        SINI)))/DJL
         B(NOINT,NOCOL3) = SML
         NOCOL4 = 3*(J-1)+M
         B(NOINT,NOCOL4) = -(SMI+SMK+SML)
  200 CONTINUE
      DOT = ZERO
      DO 220 M = 1,3
  220 DOT = DOT+C1(M)*EJI(M)/SINI
      EQVAL = PIOVR2 - ACOS(DOT)
      IF(ABS(EQVAL).LT.TOL) EQVAL=ZERO
      IF(ABS(EQVAL-PI).LT.TOL) EQVAL=PI
      RETURN
C
  800 CONTINUE
      IF (MASWRK) THEN
         WRITE(IW,9012) EJI,EJK,EJL
         WRITE (IW,9008) I,J,K,L,DOTI
      END IF
      IF(EXETYP.EQ.CHECK) RETURN
      CALL ABRT
      STOP
C
  810 CONTINUE
      IF (MASWRK) THEN
         WRITE(IW,9012) EJI,EJK,EJL
         WRITE (IW,9010) I,J,K,L,DET,SINI,SINT
      END IF
      IF(EXETYP.EQ.CHECK) RETURN
      CALL ABRT
      STOP
 9008 FORMAT(' I,J,K,L = ',4I4,5X,'DOTI =',F16.12,/' ROUNDOFF '
     *      ,'ERROR IN OPLA - PROGRAM ABORT IF EXETYP NOT CHECK')
 9010 FORMAT(' I,J,K,L = ',4I4,5X,'DET =',F16.12,5X,'SINI =',F16.12
     *      ,5X,'SINT =',F16.12,/' ROUNDOFF '
     *      ,'ERROR IN OPLA - PROGRAM ABORT IF EXETYP NOT CHECK')
 9012 FORMAT(/' EJI =',3F16.12,/' EJK =',3F16.12,/' EJL =',3F16.12)
      END
C*MODULE ZMATRX  *DECK PRZDAT
      SUBROUTINE PRZDAT(NCVAL,IZMAT,ZMAT,GRAD,IPRSW)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION IZMAT(*),ZMAT(*),GRAD(*)
      DIMENSION RTYPE(8)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /CONNEC/ IZ(MXATM,4),IATMOR(MXATM),NATM,NRATM,NDATM
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),A(MXATM),B(MXATM),BFLAB(MXAO)
C
      PARAMETER (ONE=1.0D+00, PIDEG=180.0D+00, ANGSTR=0.52917724924D+00)
C
      DATA RTYPE/8HSTRETCH ,8HBEND    ,8HTORSION ,8HPLA.BEND,
     +           8HLIN.BEND,8HDIHED. 5,8HGHOSTORS,8HSYM.CRD./
C
C     ----- PRINT THE Z-MATRIX AND IT'S GRADIENT -----
C
C           IPRSW = 0  PRINT JUST THE Z MATRIX
C                 = 1  PRINT THE Z MATRIX AND IT'S GRADIENT
C                 = 2  PRINT FULLY SUBSTITUTED Z-MATRIX
C
      IF ((IPRSW .EQ. 0) .AND. MASWRK) WRITE (IW,9008)
      IF ((IPRSW .EQ. 1) .AND. MASWRK) WRITE (IW,9028)
      IF ((IPRSW .EQ. 2) .AND. MASWRK) WRITE (IW,9038)
      RTOD = PIDEG/ACOS(-ONE)
C
      IADD = 1
      NDEGEN=0
      IF(IPRSW .EQ. 2) GO TO 300
C
      I = 0
   90 CONTINUE
      I = I + 1
      IF (I .GT. NCVAL) GO TO 210
      ITYPE = IZMAT(IADD)
C
C     ----- HANDLE SYMMETRY COORDINATES (WHICH HAVE NEGATIVE ITYPES)
C     ----- AS SEPARATE CASES
C
      IF (ITYPE .EQ. -1) THEN
C
C        ----- BOND STRETCH SYMMETRY COORDINATE -----
C
C        IPRSW = 1 IS THE ONLY CASE
C
         ZVAL=ZMAT(I)*ANGSTR
         IF (MASWRK) WRITE(IW,9109) I,RTYPE(8),ZVAL,GRAD(I)
         IADD = IADD + 1
         GO TO 200
C
      ELSE IF (ITYPE .EQ. -2) THEN
C
C        ----- ANGULAR SYMMETRY COORDINATE -----
C
C        IPRSW = 1 IS THE ONLY CASE
C
         ZVAL=ZMAT(I)*RTOD
         IF (MASWRK) WRITE(IW,9109) I,RTYPE(8),ZVAL,GRAD(I)
         IADD = IADD + 1
         GO TO 200
C
      ELSE
         GO TO (100,120,140,160,180,190,195),ITYPE
      END IF
C
C     ----- BOND STRETCH -----
C
  100 CONTINUE
      ZVAL=ZMAT(I)*ANGSTR
      IF (MASWRK) THEN
      IF (IPRSW .EQ. 0) WRITE (IW,9048) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),ZMAT(I),ZVAL
      IF (IPRSW .EQ. 1) WRITE (IW,9048) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),ZVAL,GRAD(I)
      END IF
      IADD = IADD + 3
      GO TO 200
C
C     ----- ANGLE BEND -----
C
  120 CONTINUE
      ZVAL=ZMAT(I)*RTOD
      IF (MASWRK) THEN
      IF (IPRSW .EQ. 0) WRITE (IW,9068) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),ZMAT(I),ZVAL
      IF (IPRSW .EQ. 1) WRITE (IW,9068) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),ZVAL,GRAD(I)
      END IF
      IADD = IADD + 4
      GO TO 200
C
C     ----- TORSION -----
C     ----- DIHEDRAL ANGLE BETWEEN TWO PLANES WITH TWO ATOMS IN COMMON
C
  140 CONTINUE
      ZVAL=ZMAT(I)*RTOD
      IF (MASWRK) THEN
      IF (IPRSW .EQ. 0) WRITE (IW,9088) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),IZMAT(IADD+4),ZMAT(I),ZVAL
      IF (IPRSW .EQ. 1) WRITE (IW,9088) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),IZMAT(IADD+4),ZVAL,GRAD(I)
      END IF
      IADD = IADD + 5
      GO TO 200
C
C     ----- OUT OF PLANE ANGLE BEND -----
C
  160 CONTINUE
      ZVAL=ZMAT(I)*RTOD
      IF (MASWRK) THEN
      IF (IPRSW .EQ. 0) WRITE (IW,9088) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),IZMAT(IADD+4),ZMAT(I),ZVAL
      IF (IPRSW .EQ. 1) WRITE (IW,9088) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),IZMAT(IADD+4),ZVAL,GRAD(I)
      END IF
      IADD = IADD + 5
      GO TO 200
C
C     ----- LINEAR ANGLE BEND -----
C
  180 CONTINUE
      ZVAL=ZMAT(I)*RTOD
      IF (MASWRK) THEN
      IF (IPRSW .EQ. 0) WRITE (IW,9068) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),ZMAT(I),ZVAL
      IF (IPRSW .EQ. 1) WRITE (IW,9068) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),ZVAL,GRAD(I)
      END IF
      NDEGEN=NDEGEN+1
      IF(NDEGEN.EQ.2) IADD=IADD+4
      IF(NDEGEN.EQ.2) NDEGEN=0
      GO TO 200
C
C     ----- DIHEDRAL ANGLE BETWEEN TWO PLANES WITH ONE ATOM IN COMMON
C
  190 CONTINUE
      ZVAL=ZMAT(I)*RTOD
      IF (MASWRK) THEN
      IF (IPRSW .EQ. 0) WRITE (IW,9098) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),IZMAT(IADD+4),IZMAT(IADD+5),
     *     ZMAT(I),ZVAL
      IF (IPRSW .EQ. 1) WRITE (IW,9098) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),IZMAT(IADD+4),IZMAT(IADD+5),
     *     ZVAL,GRAD(I)
      END IF
      IADD=IADD+6
      GO TO 200
C
C     ----- "GHOST" TORSION -----
C     ----- DIHEDRAL ANGLE BETWEEN TWO PLANES WITH ONE BOND IN COMMON
C
  195 CONTINUE
      ZVAL=ZMAT(I)*RTOD
      IF (MASWRK) THEN
      IF (IPRSW .EQ. 0) WRITE (IW,9108) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),IZMAT(IADD+4),IZMAT(IADD+5),
     *     IZMAT(IADD+6),ZMAT(I),ZVAL
      IF (IPRSW .EQ. 1) WRITE (IW,9108) I,RTYPE(ITYPE),IZMAT(IADD+1),
     *     IZMAT(IADD+2),IZMAT(IADD+3),IZMAT(IADD+4),IZMAT(IADD+5),
     *     IZMAT(IADD+6),ZVAL,GRAD(I)
      END IF
      IADD = IADD + 7
      GO TO 200
C
  200 CONTINUE
      GO TO 90
C
  210 CONTINUE
      RETURN
C
  300 CONTINUE
C
C PRINT FULLY SUBSTITUTED Z-MATRIX
C
      IF (.NOT.MASWRK) RETURN
C
      WRITE(IW,9110) A(IATMOR(1))
      IF(NATM.EQ.1) RETURN
      WRITE(IW,9110) A(IATMOR(2)),IZ(2,1),ZMAT(1)*ANGSTR
      IF(NATM.EQ.2) RETURN
      WRITE(IW,9110) A(IATMOR(3)),IZ(3,1),ZMAT(2)*ANGSTR,
     *                            IZ(3,2),ZMAT(3)*RTOD
      IF(NATM.EQ.3) RETURN
C
      IADD=4
      DO 400 I = 4,NATM
         WRITE(IW,9110) A(IATMOR(I)),
     *                  IZ(I,1),ZMAT(IADD)*ANGSTR,
     *                  IZ(I,2),ZMAT(IADD+1)*RTOD,
     *                  IZ(I,3),ZMAT(IADD+2)*RTOD,IZ(I,4)
         IADD = IADD + 3
  400 CONTINUE
      RETURN
C
 9008 FORMAT(//21X,20(1H-)/21X,20HINTERNAL COORDINATES/
     *         21X,20(1H-)//17X,'- - ATOMS - -',
     *         9X,'COORDINATE',6X,'COORDINATE'/
     *         1X,3HNO.,3X,4HTYPE,4X,1HI,
     *         2X,1HJ,2X,1HK,2X,1HL,2X,1HM,2X,1HN,8X,10H(BOHR,RAD),
     *         7X,9H(ANG,DEG)/1X,64(1H-))
 9028 FORMAT(//21X,20(1H-)/21X,20HINTERNAL COORDINATES/
     *         21X,20(1H-)//17X,'- - ATOMS - -',
     *         9X,'COORDINATE',8X,'GRADIENT'/
     *         1X,3HNO.,3X,4HTYPE,4X,1HI,
     *         2X,1HJ,2X,1HK,2X,1HL,2X,1HM,2X,1HN,9X,9H(ANG,DEG),
     *         5X,11H(H/B,H/RAD)/1X,64(1H-))
 9038 FORMAT(/1X,'THE CURRENT FULLY SUBSTITUTED Z-MATRIX IS')
 9048 FORMAT(1X,I3,1X,A8,2I3,9X,5X,2F16.7)
 9068 FORMAT(1X,I3,1X,A8,3I3,6X,5X,2F16.7)
 9088 FORMAT(1X,I3,1X,A8,4I3,3X,5X,2F16.7)
 9098 FORMAT(1X,I3,1X,A8,5I3,   5X,2F16.7)
 9108 FORMAT(1X,I3,1X,A8,6I3,   2X,2F16.7)
 9109 FORMAT(1X,I3,1X,A8,      20X,2F16.7)
 9110 FORMAT(1X,A4,1X,I3,F12.7,I3,F14.7,I3,F14.7,I3)
      END
C*MODULE ZMATRX  *DECK PZANDG
      SUBROUTINE PZANDG(GRAD,IPRSW)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL LINEAR,SYM
C
      DIMENSION GRAD(*)
C
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SYMINT/ SYM
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
C
C     ----- PRINT THE CURRENT Z MATRIX, AND ITS GRADIENT IF IPRSW=1-----
C
      IF (SYM  .AND.  IPRSW.EQ.1) THEN
         NCVAL = NVAR
         INZMAT = 47
         IIZMAT = 48
      ELSE
         NCVAL = NZVAR
         INZMAT = 39
         IIZMAT = 40
      END IF
C
      CALL VALFM(LOADFM)
      I10  = LOADFM + 1
      I20  = I10    + NCVAL
      LAST = I20    + NZMAT
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
C     ----- GET IZMAT AND ZMAT FROM DIRECT ACCESS -----
C     ----- PRINT THE Z MATRIX (AND GRAD) -----
C
      CALL DAREAD(IDAF,IODA,X(I10),NCVAL,INZMAT,0)
      CALL DAREAD(IDAF,IODA,X(I20),NZMAT,IIZMAT,1)
      CALL PRZDAT(NCVAL,X(I20),X(I10),GRAD,IPRSW)
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE ZMATRX  *DECK TORS
      SUBROUTINE TORS(EQVAL,NOINT,I,J,K,L,C,B,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
C     -----THIS ROUTINE COMPUTES THE B MATRIX ELEMENTS FOR THE
C          TORSION AS DEFINED BY WILSON.
C          SEE PP. 60-61 OF "MOLECULAR VIBRATIONS" BY E.B. WILSON,
C          J.C. DECIUS AND P.C. CROSS, MCGRAW-HILL, NY, 1955.
C          NOTE THAT THE SIGN OF THE SECOND TERM ON THE RIGHT HAND
C          SIDE IN EQUATION (22) IS INCORRECT AS PRINTED ON P61 OF
C          WDC -----
C
C     -----I AND L ARE THE END ATOMS.  I IS NEARER THE OBSERVER.
C          J AND K ARE THE CENTRAL ATOM. J IS NEARER THE OBSERVER -----
C
      DIMENSION C(3,1),B(NDIM,1)
      DIMENSION RIJ(3),RJK(3),RKL(3),EIJ(3),EJK(3),EKL(3),CR1(3),CR2(3)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TOL=1.0D-07)
      PARAMETER (TOLRD=1.0001D+00)
C
C
      PI = ACOS(-ONE)
      DIJSQ = ZERO
      DJKSQ = ZERO
      DKLSQ=ZERO
      DO 120 M = 1,3
         RIJ(M) = C(M,J)-C(M,I)
         DIJSQ = DIJSQ+RIJ(M)*RIJ(M)
         RJK(M) = C(M,K)-C(M,J)
         DJKSQ = DJKSQ+RJK(M)*RJK(M)
         RKL(M) = C(M,L)-C(M,K)
         DKLSQ = DKLSQ+RKL(M)*RKL(M)
  120 CONTINUE
      DIJ = SQRT(DIJSQ)
      DJK = SQRT(DJKSQ)
      DKL = SQRT(DKLSQ)
      DO 180 M = 1,3
         EIJ(M) = RIJ(M)/DIJ
         EJK(M) = RJK(M)/DJK
         EKL(M) = RKL(M)/DKL
  180 CONTINUE
C
      CR1(1) = EIJ(2)*EJK(3)-EIJ(3)*EJK(2)
      CR1(2) = EIJ(3)*EJK(1)-EIJ(1)*EJK(3)
      CR1(3) = EIJ(1)*EJK(2)-EIJ(2)*EJK(1)
C
      CR2(1) = EJK(2)*EKL(3)-EJK(3)*EKL(2)
      CR2(2) = EJK(3)*EKL(1)-EJK(1)*EKL(3)
      CR2(3) = EJK(1)*EKL(2)-EJK(2)*EKL(1)
C
      DOTPJ = ZERO
      DOTPK = ZERO
      DO 220 M = 1,3
         DOTPJ = DOTPJ-EIJ(M)*EJK(M)
         DOTPK = DOTPK-EJK(M)*EKL(M)
  220 CONTINUE
      IF (ONE.LE.ABS(DOTPJ)) GO TO 320
      IF (ONE.LE.ABS(DOTPK)) GO TO 320
      SINPJ = SQRT(ONE-DOTPJ*DOTPJ)
      SINPK = SQRT(ONE-DOTPK*DOTPK)
      DO 280 M = 1,3
         SMI = -CR1(M)/(DIJ*SINPJ*SINPJ)
         NOCOL1 = 3*(I-1)+M
         B(NOINT,NOCOL1) = SMI
         F1 = (CR1(M)*(DJK-DIJ*DOTPJ))/(DJK*DIJ*SINPJ*SINPJ)
         F2 = (DOTPK*CR2(M))/(DJK*SINPK*SINPK)
         SMJ = F1-F2
         NOCOL2 = 3*(J-1)+M
         B(NOINT,NOCOL2) = SMJ
         SML = CR2(M)/(DKL*SINPK*SINPK)
         NOCOL3 = 3*(L-1)+M
         B(NOINT,NOCOL3) = SML
         NOCOL4 = 3*(K-1)+M
         B(NOINT,NOCOL4) = -(SMI+SMJ+SML)
  280 CONTINUE
      DOT = ZERO
      DO 300 M = 1,3
  300 DOT = DOT+(CR1(M)*CR2(M))/(SINPJ*SINPK)
C
      IF(ABS(DOT).LE.TOLRD) THEN
         IF (ABS(DOT) .GT. ONE) DOT = SIGN(ONE,DOT)
         EQVAL = ACOS(DOT)
         IF(ABS(EQVAL).LT.TOL) EQVAL=ZERO
         IF(ABS(EQVAL-PI).LT.TOL) EQVAL=PI
         SENSE = DDOT(3,CR2,1,RIJ,1)
         IF (SENSE.LT.ZERO .AND.  ABS(SENSE).GT.TOL) EQVAL = -EQVAL
         RETURN
      END IF
C
  320 CONTINUE
      IF (MASWRK) WRITE (IW,9008)
      CALL ABRT
      STOP
 9008 FORMAT(1X,'ARITHMETIC ROUNDOFF ERROR IN TORS - STOP')
      END
C*MODULE ZMATRX  *DECK ZMATIN
      SUBROUTINE ZMATIN(CARTSN)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION IDUMMY(2)
      PARAMETER (MXATM=2000)
      LOGICAL CARTSN
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
C     ----- GROW FAST MEMORY FOR $ZMAT READING -----
C
      CALL GTBOND(IDUMMY,1,NBONDS,-1)
      NPAIR =    NBONDS +  20
      NANG  =  5*NBONDS + 100
      NTOR  = 20*NBONDS + 800
      NNVDW = 2*NAT
      NELS = 30*MAX(3*NAT-6,1)
      NSIJ = (NPAIR+NANG+NTOR) * 3 * NAT
      NIJS = 10*3*NAT
C
C     DLC : ADDED NPAIR : NUMBER OF BONDS
C                 NANG  : NUMBER OF ANGLES
C                 NTOR  : NUMBER OF TORSIONAL ANGLES
C                 NNVDW : NUMBER OF VAN DER WAALS PAIRS
C                 LNXZMAT : TEM. STORAGE FOR EXTRA INTERNAL COORDINATES.
C                 LNRZMAT : TEM. STORAGE FOR REMOVED INTERNAL COORDS
C                 LNFZMAT : TEM. STORAGE FOR CONSTRAINED INTERNAL COORDS
C                 LNTZMAT : TEM. STORAGE FOR TEM. INTERNAL COORDINATES.
C                           IT IS NEEDED TO REORDER THE INTERNAL COORDS.
C
      CALL VALFM(LOADFM)
      LIZMAT  = LOADFM  + 1
      LSIJ    = LIZMAT  + NELS
      LIJS    = LSIJ    + NSIJ
      LNVDW   = LIJS    + 2*NIJS
      LNXZMAT = LNVDW   + NNVDW
      LNFZMAT = LNXZMAT + NELS
      LNTZMAT = LNFZMAT + NELS
      LNRZMAT = LNTZMAT + NELS
      LNPAIR  = LNRZMAT + NELS
      LNANG   = LNPAIR  + 2*NPAIR
      LNTOR   = LNANG   + 3*NANG
      LAST    = LNTOR   + 4*NTOR
C
      NEED = LAST - LOADFM
      CALL GETFM(NEED)
C
      CALL ZMAT2(CARTSN,X(LSIJ),X(LIJS),X(LIZMAT),X(LNVDW),X(LNXZMAT),
     *     X(LNFZMAT),X(LNTZMAT),X(LNRZMAT),X(LNPAIR),X(LNANG),X(LNTOR),
     *     NELS,NIJS,NSIJ,NPAIR,NANG,NTOR,NNVDW)
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE ZMATRX  *DECK ZMAT2
      SUBROUTINE ZMAT2(CARTSN,SIJ,IJS,IZMAT,NONVDW,IXZMAT,IFZMAT,
     *                 ITZMAT,IRZMAT,IPAIR,IANG,ITOR,NELS,
     *                 NIJS,NSIJ,NPAIR,NANG,NTOR,NNVDW)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION IZMAT(NELS),IJS(2,NIJS),SIJ(NSIJ)
      DIMENSION IXZMAT(NELS),IFZMAT(NELS),ITZMAT(NELS),NONVDW(NNVDW),
     *          IRZMAT(NELS),IPAIR(2,NPAIR),IANG(3,NANG),ITOR(4,NTOR)
C
      LOGICAL LINEAR,DECOMP,NVAROK,GOPARR,DSKWRK,MASWRK
      LOGICAL SYM,DLC,AUTO,CARTSN
      LOGICAL MMONLY,QMMM,IMOMM,SIMOMM
C
      PARAMETER (MXATM=2000, MAXCST=50, MAXLNK=50)
C
      COMMON /DLCFRZ/ FVALUE(MAXCST),ITABLE(MAXCST),IFTYPE(MAXCST),
     *                NCONST
      COMMON /INFOA / NAT,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 /QMMM1 / IMOMM,SIMOMM,NAPAIR,NSEQ
      COMMON /QMMM3 / LINKGE(3*MAXLNK),BLQM(MAXLNK),BLMM(MAXLNK)
      COMMON /SYMINT/ SYM
      COMMON /TINOPT/ MMONLY,QMMM
      COMMON /VBDCMP/ DECOMP,NROW,NDEG
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
      COMMON /ZMTALT/ NZMAT2,NZVAR2,NVAR2,NZMTRD,ICOORD
C
      PARAMETER (ZERO=0.0D+00)
C
      PARAMETER (NNAM=10)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA SZMAT/8HZMAT    /
      DATA QNAM /8HIZMAT   ,8HSIJ     ,8HIJS     ,8HDLC     ,8HAUTO    ,
     *           8HNONVDW  ,8HIXZMAT  ,8HIFZMAT  ,8HIRZMAT  ,8HFVALUE  /
C
      DATA KQNAM/-1,-3,-1,0,0,1,1,1,1,-3/
C
C     ----- THIS ROUTINE READS THE CODED Z-MATRIX -----
C     IN COMMON /ZMAT/, NZVAR IS THE NUMBER OF INPUT
C     COORDINATES, NVAR IS THE TOTAL NUMBER OF COORDINATES,
C     AND NZMAT IS THE TOTAL NUMBER OF CODE INTEGERS INPUT.
C
      KQNAM(1)  = 10*NELS   + 1
      KQNAM(2)  = 10*NSIJ   + 3
      KQNAM(3)  = 10*2*NIJS + 1
      KQNAM(6)  = 10*NNVDW  + 1
      KQNAM(7)  = 10*NELS   + 1
      KQNAM(8)  = 10*NELS   + 1
      KQNAM(9)  = 10*NELS   + 1
      KQNAM(10) = 10*MAXCST + 3
C
      DO I=1,NELS
         IZMAT(I) = 0
         IXZMAT(I) = 0
         IFZMAT(I) = 0
         IRZMAT(I) = 0
      ENDDO
      DO I=1,MAXCST
         FVALUE(I) = ZERO
      ENDDO
      DO I = 1,NIJS
         IJS(1,I) = 0
         IJS(2,I) = 0
      ENDDO
      DO I = 1,NSIJ
         SIJ(I) = ZERO
      ENDDO
      DLC = .FALSE.
      AUTO= .FALSE.
      JRET=0
      DO I=1,NNVDW
         NONVDW(I) = 0
      ENDDO
C
      CALL NAMEIO(IR,JRET,SZMAT,NNAM,QNAM,KQNAM,
     *            IZMAT,SIJ,IJS,DLC,AUTO,NONVDW,
     *            IXZMAT,IFZMAT,IRZMAT,FVALUE,
     *            0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
C
      IF ((DLC .OR. AUTO)  .AND.  LINEAR) THEN
         IF(MASWRK) THEN
            WRITE(IW,*) 'DLC=.TRUE. AND AUTO=.TRUE. ARE NOT ',
     *                  'COMPATIBLE WITH LINEAR MOLECULES.'
            WRITE(IW,*) 'PLEASE USE LINEAR BEND COORDINATES IN $ZMAT.'
         END IF
         CALL ABRT
      END IF
C
C ---- AUTOMATICALLY GENERATE A SET OF INTERNAL COORDIANTES
C
      IF (DLC .AND. AUTO) THEN
C
C    LET'S ADD DEFAULT GEOMETRY CONSTRAINTS OF IMOMM
C
         IF (QMMM) THEN
            IF ((IMOMM .AND. (ICOORD.EQ.5)) .OR.
     *          (IMOMM .AND. (ICOORD.EQ.-1)) ) THEN
               ILAST=0
               DO I=1,MAXCST
                  IF (FVALUE(I).NE.ZERO) ILAST=ILAST+1
               ENDDO
               IZLAST=0
               DO 450 I=1,ILAST
                  IF (IFZMAT(IZLAST).EQ.1) THEN
                     IZLAST=IZLAST+3
                  ELSE IF (IFZMAT(IZLAST).EQ.2) THEN
                     IZLAST=IZLAST+4
                  ELSE
                     IZLAST=IZLAST+5
                  ENDIF
 450           CONTINUE
               DO I=1,NAPAIR
                  FVALUE(I+ILAST)=BLQM(I)
                  IFZMAT(IZLAST+I*3-2)=1
                  IFZMAT(IZLAST+I*3-1)=LINKGE(I*3-2)
                  IFZMAT(IZLAST+I*3  )=LINKGE(I*3-1)
               ENDDO
            ENDIF
         ENDIF
         CALL AUTOIZ(IZMAT,NONVDW,IXZMAT,IFZMAT,ITZMAT,IRZMAT,
     *               IPAIR,IANG,ITOR,NELS,NNVDW,NZVAR,NPAIR,NANG,NTOR)
         GOTO 200
C
C          LET NON-AUTO RUNS HAVE ACCESS TO FROZEN COORDINATE CAPABILITY
C
      ELSE
         NXCONST=0
C
C    1. IS INTERNAL COORDINATE OF IFZMAT IN IZMAT ?
C    2. FIND OUT THE NUMBER IN IZMAT OF THE EACH INTERNAL COORDINATE
C       OF IFZMAT.
C    3. GET NCONST
C
         CALL CHKFZT(IFZMAT,IZMAT,NELS,NZVAR)
C
C -- PUT IFZMAT (CONSTRAINT INTERNAL COORDINATE) INTO
C -- IZMAT ARRAY, IF IT IS NOT DEFINED.
C
         IF (NCONST.GT.0) THEN
            NZMAT=0
            DO I=1,NZVAR
               ITYPE = IZMAT(NZMAT+1)
               IF(ITYPE.EQ.1) NZMAT=NZMAT+3
               IF(ITYPE.EQ.2) NZMAT=NZMAT+4
               IF(ITYPE.EQ.3) NZMAT=NZMAT+5
               IF(ITYPE.EQ.4) NZMAT=NZMAT+5
               IF(ITYPE.EQ.5) NZMAT=NZMAT+4
               IF(ITYPE.EQ.6) NZMAT=NZMAT+6
               IF(ITYPE.EQ.7) NZMAT=NZMAT+7
            ENDDO
            CALL ICOPY(NZMAT,IZMAT,1,ITZMAT,1)
            IZ = NZMAT+1
            NIZ = 1
            DO 840 I=1,NCONST
              IF (IFZMAT(NIZ) .EQ. 1) THEN
                 IF (ITABLE(I) .EQ. -1 ) THEN
                   ITZMAT(IZ  ) = IFZMAT(NIZ)
                   ITZMAT(IZ+1) = IFZMAT(NIZ+1)
                   ITZMAT(IZ+2) = IFZMAT(NIZ+2)
                   IZ = IZ + 3
                   NIZ = NIZ + 3
                   NZVAR = NZVAR + 1
                   NZMAT = NZMAT + 3
                   ITABLE(I) = NZVAR
                   NXCONST = NXCONST + 1
                 ELSE
                   NIZ = NIZ + 3
                 END IF
              ELSEIF (IFZMAT(NIZ) .EQ. 2) THEN
                 IF (ITABLE(I) .EQ. -1 ) THEN
                   ITZMAT(IZ  ) = IFZMAT(NIZ)
                   ITZMAT(IZ+1) = IFZMAT(NIZ+1)
                   ITZMAT(IZ+2) = IFZMAT(NIZ+2)
                   ITZMAT(IZ+3) = IFZMAT(NIZ+3)
                   IZ = IZ + 4
                   NIZ = NIZ + 4
                   NZVAR = NZVAR + 1
                   NZMAT = NZMAT + 4
                   ITABLE(I) = NZVAR
                   NXCONST = NXCONST + 1
                 ELSE
                   NIZ = NIZ + 4
                 END IF
              ELSEIF (IFZMAT(NIZ) .EQ. 3) THEN
                 IF (ITABLE(I) .EQ. -1 ) THEN
                   ITZMAT(IZ  ) = IFZMAT(NIZ)
                   ITZMAT(IZ+1) = IFZMAT(NIZ+1)
                   ITZMAT(IZ+2) = IFZMAT(NIZ+2)
                   ITZMAT(IZ+3) = IFZMAT(NIZ+3)
                   ITZMAT(IZ+4) = IFZMAT(NIZ+4)
                   IZ = IZ + 5
                   NIZ = NIZ + 5
                   NZVAR = NZVAR + 1
                   NZMAT = NZMAT + 5
                   ITABLE(I) = NZVAR
                   NXCONST = NXCONST + 1
                 ELSE
                   NIZ = NIZ + 5
                 END IF
              ENDIF
 840        CONTINUE
C
C     LET'S REORDER THE IZMAT IN SUCH A WAY THAT CONSTRAINT
C     COORDINATES COME FIRST.
C
            CALL RRD(IZMAT,ITZMAT,NELS,NZVAR)
         ENDIF
      ENDIF
C
      NZMTRD = JRET
      IF(JRET.EQ.0) GO TO 200
      IF(JRET.EQ.2) THEN
         IF (MASWRK) WRITE (IW,9008)
         CALL ABRT
      END IF
C
C         INPUT OF INTERNALS REQUIRED IF COORD=UNIQUE, HINT, OR CART,
C         BUT FOR COORD=ZMT, ZMTMPC WE CAN USE $DATA INFO INSTEAD
C
      IF(JRET.EQ.1  .AND.  ICOORD.LE.1) THEN
         IF (MASWRK) WRITE(IW,9008)
         CALL ABRT
      ELSE
         CALL FILLZM(IZMAT)
      END IF
C
C     ----- SCAN AND ANALYZE IZMAT -----
C     THE VALID TYPES ARE
C        1= STRETCH I-J
C        2= BEND I-J-K
C        3= TORSION I-J-K-L
C           (BETWEEN TWO PLANES HAVING TWO ATOMS IN COMMON)
C        4= OUT OF PLANE BEND I-J-K-L
C        5= LINEAR BEND I-J-K
C        6= 5 ATOM DIHEDRAL ANGLE I-J-K-L-M
C           (BETWEEN TWO PLANES HAVING ONE ATOM IN COMMON)
C        7= 6 ATOM GHOST TORSION ANGLE I1-I2-J-K-L1-L2
C
  200 CONTINUE
      IF (MASWRK) WRITE(IW,9048)
      IF (MASWRK) WRITE(IP,8000)
      NVAR=0
      IZ = 1
      DO 230 I=1,NZVAR
         ITYPE=IZMAT(IZ)
         NCHK=0
         IF(ITYPE.EQ.1) NCHK=2
         IF(ITYPE.EQ.2) NCHK=3
         IF(ITYPE.EQ.3) NCHK=4
         IF(ITYPE.EQ.4) NCHK=4
         IF(ITYPE.EQ.5) NCHK=3
         IF(ITYPE.EQ.6) NCHK=5
         IF(ITYPE.EQ.7) NCHK=6
         IF(NCHK.EQ.0) THEN
            IF (MASWRK) WRITE(IW,9068) I,IZ,IZMAT(IZ)
            CALL ABRT
         END IF
C
         IF (MASWRK) WRITE(IW,9058) I,ITYPE,(IZMAT(IZ+J),J=1,NCHK)
         IF (MASWRK) WRITE(IP,8010) ITYPE,(IZMAT(IZ+J),J=1,NCHK)
         DO 210 J=1,NCHK
            IAT=IZMAT(IZ+J)
            IF(IAT.LE.0  .OR.  IAT.GT.NAT) THEN
               IF (MASWRK) WRITE(IW,9088) I,IZ,IAT,J,NAT
               CALL ABRT
            END IF
  210    CONTINUE
C
         IZ=IZ+NCHK+1
         NVAR=NVAR+1
         IF(ITYPE.EQ.5) NVAR=NVAR+1
  230 CONTINUE
      NZMAT = IZ - 1
      NZVAR = NVAR
C
      IF (MASWRK) WRITE(IP,8020)
C
C --- IF NOT DOING VIBRATIONAL ANALYSIS, MUST SPECIFY EXACTLY 3N-NDEG
C --- INTERNAL COORDINATES.  OTHERWISE, MORE THAN 3N-NDEG IS OK.
C
      NDEG = 6
      IF (LINEAR) NDEG = 5
      N3NM6 = 3*NAT-NDEG
      NVAROK = NVAR .EQ. (3*NAT-NDEG)
C
C     ----- SAVE IZMAT ON DIRECT ACCESS -----
C     DLC : SINCE WE GENERATED REDUNDANT INTERNAL COORDINATES.
C
      CALL DAWRIT(IDAF,IODA,IZMAT,NZMAT,40,1)
C
C     ----- READ IN SYMMETRIC COORDINATES MATRIX -----
C
      CALL ZMTSMM(SIJ,IJS,IZMAT,NZVAR,N3NM6,NIJS,NSIJ,NZMAT,NSYMC,
     +            SYM,DLC)
C
C     ----- IF SYMMETRY COORDINATES ARE DEFINED, DISABLE
C           VIBRATIONAL DECOMPOSITION ------
C
      IF (NSYMC.GT.0 .AND. DECOMP) THEN
         DECOMP = .FALSE.
         IF (MASWRK) WRITE(IW,9208)
      END IF
C
      IF (CARTSN) NVAROK=.TRUE.
      IF (DECOMP .OR. SYM) NVAROK = NVAR .GE. (3*NAT-NDEG)
      IF (.NOT.NVAROK) THEN
         IF (MASWRK) WRITE(IW,9148) NDEG,NAT,NVAR
         CALL ABRT
      END IF
C
C     ----- WE CAN ONLY OPTIMIZE IN 3N-6(5) COORDINATES, SO UNLESS WE
C           ARE DOING A DECOMPOSITION NVAR = 3N-6(5) -----
C
      IF (.NOT. DECOMP) NVAR = N3NM6
      NROW = NVAR
C
C     ----- SAVE IZMAT ON DIRECT ACCESS -----
C
      CALL DAWRIT(IDAF,IODA,IZMAT,NZMAT,40,1)
C
      NVAR2  = NVAR
      NZMAT2 = NZMAT
      RETURN
C
 8000 FORMAT(1X,'$ZMAT   IZMAT(1)=')
 8010 FORMAT(6X,7(I4,','))
 8020 FORMAT(1X,'$END')
C
 9008 FORMAT(1X,'ERROR IN $ZMAT INPUT - ERROR')
 9048 FORMAT(/3X,'--- ENCODED Z MATRIX ---'/
     *        1X,'COORD',2X,'TYPE',3X,'I',3X,'J',3X,'K',3X,'L',3X,'M',
     *        3X,'N')
C     DLC : INITIAL NUMBER OF REDUANDANT INTERNAL COORDINATE
C           CAN BE LARGE. SO I3 -> I4
 9058 FORMAT(1X,I4,3X,I4,1X,6I4)
 9068 FORMAT(1X,'INTERNAL COORD.',I3,' AT SUBSCRIPT',I4,
     *     ' IS OF ILLEGAL TYPE',I5)
 9088 FORMAT(1X,'INTERNAL COORD.',I3,' AT SUBSCRIPT',I4,
     *    ' IS FOLLOWED BY ILLEGAL ATOM',I4,' AT OFFSET',I4/
     *    ' HIGHEST LEGAL ATOM NUMBER IS',I4)
 9148 FORMAT(1X,'YOU MUST INPUT 3N-',I1,' INTERNAL COORDS WHERE N ='
     *      ,I3/' YOU HAVE PROVIDED ',I6)
 9208 FORMAT(/1X,'*** SYMMETRY COORDINATES ARE DEFINED, SO ',
     *          'VIBRATIONAL DECOMPOSITION IS DISABLED ***')
      END
C*MODULE ZMATRX  *DECK ZMTSMM
      SUBROUTINE ZMTSMM(SIJ,IJS,IZMAT,NZVAR,N3NM6,NIJS,NSIJ,NZMAT,
     *                  NSYMC,SYM,DLC)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SYM,DLC
C
      COMMON /FMCOM / X(1)
C
      DIMENSION IJS(2,NIJS), SIJ(NSIJ), IZMAT(*)
C
C     ----- GROW FAST MEMORY -----
C
      NDUM = NZVAR*NZVAR
      CALL VALFM(LOADFM)
      LSYM   = LOADFM + 1
      LIZMT2 = LSYM   + NDUM
      LAST   = LIZMT2 + NZMAT
      NEED   = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      CALL ZMTSYM(X(LSYM),SIJ,IJS,IZMAT,X(LIZMT2),NZVAR,
     *            N3NM6,NZMAT,NIJS,NSIJ,NSYMC,SYM,DLC)
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE ZMATRX  *DECK ZMTSYM
      SUBROUTINE ZMTSYM(S,SIJ,IJS,IZMAT,IZMAT2,NZVAR,
     *                  N3NM6,NZMAT,NIJS,NSIJ,NSYMC,SYM,DLC)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION IJS(2,NIJS),SIJ(NSIJ),S(NZVAR*NZVAR),IZMAT(*),
     *          IZMAT2(NZMAT)
      PARAMETER (LBUF = 5)
      DIMENSION BUF(LBUF),IBUF(LBUF),ISTP(8)
C
      LOGICAL SYM, DLC
      LOGICAL OUT, MASWRK, DSKWRK, GOPARR
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TOL=1.0D-08, TENM1=1.0D-01)
C
      DATA ISTP/3,4,5,5,4,6,7,1/
C
      OUT = NPRINT .NE. -5  .AND.  MASWRK
C
C     ----- NSYMC IS THE NUMBER OF SYMMETRY COORDINATES.
C
      NSYMC = 0
      NDUM = NZVAR*NZVAR
      CALL VCLR(S,1,NDUM)
      DO 20 I = 1,NZVAR-N3NM6
         II = NZVAR*(N3NM6-1 + I) + (N3NM6 + I)
         S(II) = ONE
   20 CONTINUE
C
C     ----- SET UP SYMMETRY COORDINATE TRANSFORMATION MATRIX -----
C
      IF(DLC) THEN
         SYM = .TRUE.
C     DLC : CALCULATE COEFFICIENTS OF DLCS
         CALL DLCBGD(SIJ,NSIJ)
C
         N = 0
         DO 250 J = 1,N3NM6
            DO 240 I = 1,NZVAR
               N = N+1
               IJ = NZVAR*(J-1) + I
               S(IJ) = SIJ(N)
  240       CONTINUE
  250    CONTINUE
      ELSE
         SYM = .FALSE.
         DO 370 N = 1,NIJS
            I = IJS(1,N)
            J = IJS(2,N)
            SVAL = SIJ(  N)
C  ARE ALL VALUES BETWEEN 0-NZVAR?
            IF (      I) 820,370,330
  330       IF (NZVAR-I) 820,340,340
  340       IF (      J) 820,370,350
  350       IF (NZVAR-J) 820,360,360
  360       CONTINUE
            IJ = NZVAR*(J-1) + I
            S(IJ) = SVAL
            SYM = .TRUE.
  370    CONTINUE
      END IF
C
      IF (.NOT.SYM) GO TO 600
C
C     ----- NORMALIZE THE COLUMNS OF S -----
C
      DO 470 J = 1,NZVAR
         DUM = ZERO
         DO 440 I = 1,NZVAR
            IJ = NZVAR*(J-1) + I
            DUM = DUM + S(IJ)**2
  440    CONTINUE
         IF (DUM .LT. TOL) GO TO 840
         DUM = ONE/SQRT(DUM)
         DO 460 I = 1,NZVAR
            IJ = NZVAR*(J-1) + I
            S(IJ) = S(IJ)*DUM
  460    CONTINUE
  470 CONTINUE
C
C     ----- PRINT SYMMETRY MATRIX -----
C
      IF (OUT) WRITE(IW,9997)
      IF (.NOT. OUT) GO TO 700
      DO 520 J=1,N3NM6
         NBUF=1
         N=0
         DO 510 I=1,NZVAR
            IJ=NZVAR*(J-1)+I
            DUM=S(IJ)
            IF( ABS(DUM).LT.TENM1) GO TO 510
            N=N+1
            IBUF(N)=I
            BUF(N)=DUM
            IF (N.LT.LBUF) GO TO 510
            IF (NBUF.EQ.1) WRITE(IW,9995) J,(BUF(M),IBUF(M),M=1,N)
            IF (NBUF.GT.1) WRITE(IW,9994)   (BUF(M),IBUF(M),M=1,N)
            NBUF=NBUF+1
            N=0
  510    CONTINUE
         IF (N.EQ.0) GO TO 520
         IF (NBUF.EQ.1) WRITE(IW,9995) J,(BUF(M),IBUF(M),M=1,N)
         IF (NBUF.GT.1) WRITE(IW,9994)   (BUF(M),IBUF(M),M=1,N)
  520 CONTINUE
      GO TO 700
C
C     ----- NO SYMMETRY COORDINATES GIVEN: S = I -----
C
  600 CONTINUE
      IF (NZVAR.GT.N3NM6) WRITE(IW,9996)
      DO 610 I = 1,NZVAR
         II = NZVAR*(I-1) + I
         S(II) = ONE
  610 CONTINUE
C
C     ----- GENERATE NEW IZMAT -----
C     IADD,  NDEG  = COUNTER FOR OLD IZMAT
C     IADD1, NDEG1 = COUNTER FOR NEW IZMAT2
C
  700 CONTINUE
      IADD1 = 1
      NDEG1 = 0
      DO 750 J = 1,N3NM6
         DUM = ZERO
C
C     ----- IS INTERNAL COORD. J A SYM. COORD.?
C
         DO 720 K = 1,NZVAR
            IJ = NZVAR*(J-1) + K
            DUM = DUM + ABS( S(IJ) )
            IF (ABS(S(IJ)-ONE).LT.1.0D-08) ICRD = K
  720    CONTINUE
C
         IF(ABS(DUM-ONE).LT.1.0D-08) THEN
C
C                SKIP PAST COORDINATES PRECEEDING -ICRD-
C
            NDEG = 0
            IADD = 1
            DO 730 K = 1,ICRD-1
               ITYPE = IZMAT(IADD)
               IF(ITYPE.EQ.5) THEN
                  NDEG = NDEG+1
                  IF(NDEG.EQ.1) IADD = IADD - ISTP(ITYPE)
                  IF(NDEG.EQ.2) NDEG=0
               END IF
               IADD = IADD + ISTP(ITYPE)
  730       CONTINUE
C
            ITYPE = IZMAT(IADD)
            IF(ITYPE.LE.0  .OR.  ITYPE.GT.7) THEN
               IF(MASWRK) WRITE(IW,9990) J,ITYPE,IADD
               CALL ABRT
            END IF
            IZMAT2(IADD1) = ITYPE
C
            KMAX = ISTP(ITYPE)-1
            DO 740 K = 1,KMAX
               IADD1 = IADD1 + 1
               IZMAT2(IADD1) = IZMAT(IADD+K)
  740       CONTINUE
            IF(ITYPE.EQ.5) THEN
               NDEG1=NDEG1+1
               IF(NDEG1.EQ.1) IADD1 = IADD1 - ISTP(ITYPE)
               IF(NDEG1.EQ.2) NDEG1=0
            END IF
         ELSE
C
C           ----- ALL SYMMETRY COORDINATES ARE INITIALLY DESIGNATED
C           ----- AS TYPE -1; LATER ON ALL NON-STRETCHING SYM. COORDS
C           ----- WILL BE CHANGED TO TYPE -2.
C
C
            IZMAT2(IADD1) = -1
            NSYMC = NSYMC + 1
         END IF
         IADD1 = IADD1 + 1
  750 CONTINUE
C
C     ----- IF SYMMETRY COORDINATES ARE USED, DETERMINE WHICH ARE
C     ----- STRETCHES AND WHICH ARE ANGLES.
C
      IF(NSYMC.GT.0  .AND.  .NOT.DLC) THEN
         IADD1 = 1
         NDEG1 = 0
         IJSPT = 0
         DO 780 J=1,N3NM6
            IJSPT = IJSPT+1
            ITYP1 = IZMAT2(IADD1)
C
C           ----- SKIP OVER "PRIMITIVE" INTERNAL COORDINATES
C
            IF (ITYP1.GT.0) THEN
               IADD1 = IADD1 + ISTP(ITYP1)
               IF (ITYP1.EQ.5) THEN
                  NDEG1 = NDEG1 + 1
                  IF (NDEG1.EQ.1) IADD1 = IADD1 - ISTP(ITYP1)
                  IF (NDEG1.EQ.2) NDEG1 = 0
               END IF
C
C           ----- IS THIS A BOND-STRETCH SYMMETRY COORDINATE?
C
            ELSE
C
C           ----- IPRIM IS THE "IPRIM"TH PRIMITIVE INTERNAL COORDINATE.
C           ----- JREF IS ONE OF THE SYMMETRY COORDINATES TO WHICH THE
C           ----- "IPRIM"TH PRIMITIVE INTERNAL COORDINATE CONTRIBUTES.
C
               IPRIM = IJS(1,IJSPT)
C---           JREF  = IJS(2,IJSPT)
C
C              ----- DETERMINE THE TYPE OF INTERNAL COORDINATE "IPRIM"
C
               NDEG = 0
               IADD = 1
               DO 760 K = 1,IPRIM-1
                  ITYPE = IZMAT(IADD)
                  IF(ITYPE.EQ.5) THEN
                     NDEG = NDEG + 1
                     IF(NDEG.EQ.1) IADD = IADD - ISTP(ITYPE)
                     IF(NDEG.EQ.2) NDEG=0
                  END IF
                  IADD = IADD + ISTP(ITYPE)
  760          CONTINUE
C
               ITYPE = IZMAT(IADD)
               IF (ITYPE.GT.1) IZMAT2(IADD1) = -2
C
C              ----- ADVANCE TO THE NEXT COORDINATE BY INCREMENTING
C              ----- IJSPT UNTIL IJS(2,IJSPT) CHANGES.
C
  770          CONTINUE
               IF (IJS(2,IJSPT+1).EQ.IJS(2,IJSPT)) THEN
                  IJSPT = IJSPT+1
                  GO TO 770
C                 *********
               END IF
               IADD1 = IADD1 + 1
            END IF
  780    CONTINUE
      END IF
C
C     ----- SAVE SYMMETRY MATRIX -----
C
      NDUM = NZVAR*NZVAR
C
      CALL DAWRIT(IDAF,IODA,S,NDUM,46,0)
      CALL DAWRIT(IDAF,IODA,IZMAT2,NZMAT,48,1)
      RETURN
C
C     ----- ERRORS -----
C
  820 CONTINUE
      IF(MASWRK) WRITE(IW,9999) I,J,S
      CALL ABRT
      STOP
  840 CONTINUE
      IF(MASWRK) WRITE(IW,9998) J
      CALL ABRT
      STOP
C
 9999 FORMAT(' INPUT ERROR IN -ZMTSYM- . I,J,S = ',2I5,F10.5)
 9998 FORMAT(' NORMALIZATION ERROR IN -ZMTSYM- FOR VECTOR J = ',I5)
 9997 FORMAT(/,10X,'INDEP./SYM. INTERNAL COORDINATES',/,10X,32(1H-))
 9996 FORMAT(' WARNING. NUMBER OF INTERNAL COORDINATES IS GREATER',
     * ' THAN (3N-6), BUT NO SYMMETRY COORDINATES ARE GIVEN.')
 9995 FORMAT(1X,'CRD',I3,'= ',5(F8.4,' (',I3,')'))
 9994 FORMAT(9X,5(F8.4,' (',I3,')'))
 9990 FORMAT(1X,'ERROR BUILDING NON-REDUNDANT -IZMAT-'/
     *       1X,'CHECK LINEAR COMBINATION',I5/
     *       1X,'TYPE=',I5,' ??? AT POSITION',I5,' IN YOUR INPUT IZMAT')
      END
C*MODULE ZMATRX  *DECK SQZMAT
      SUBROUTINE SQZMAT(VEC,ICOL,IROW,LDV)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION VEC(LDV*ICOL)
C
      IF (IROW .EQ. LDV) RETURN
C
      DO 100 IIP = 1,ICOL-1
         IS =     IIP*LDV + 1
         ISHDBE = IIP*IROW + 1
         CALL DCOPY(IROW,VEC(IS),1,VEC(ISHDBE),1)
  100 CONTINUE
C
      RETURN
      END
