C 14 NOV 05 - DGF - PREVS: STARTING MO ARG ADDED
C 19 SEP 05 - IA  - SYNCHRONIZE FRGINF COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  1 JUN 05 - MWS - PRATM: ENSURE SERIAL PRINTING ONLY
C 30 APR 05 - DGF - SYNCHRONISE SYMBLK COMMON BLOCK
C 14 MAR 05 - HL  - CHANGE FMASS IN FRGMSS COMMON
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - ADD ROUTINE TO INITIALIZE GAUSS-HERMITE QUADRATURES
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C  7 APR 04 - RMO - INCLUDE POLYNOMIAL FIT ROUTINE
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 14 JAN 03 - JAB - DGECO: REMOVE DOUBLE PRECISION FUNCTION
C 17 APR 02 - MWS - SYNCH UP EFMULT AND FRGINF COMMON
C 26 MAR 02 - MWS - INCLUDE LINPACK ROUTINE DGECO
C 16 NOV 01 - JMS - MOVE GAMGEN TO INT2B MODULE
C  8 OCT 01 - JMS - GAMGEN: CUG1=16 INCREASED TO 18
C 19 SEP 01 - MWS - CONVERT MXAOCI PARAMETER TO MXAO
C  1 AUG 01 - JMS - GAMGEN: REVISED FOR NEW ROTATED AXIS INT/GRAD CODES
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 13 JUN 01 - DGF - SEREP: BAIL OPTION ADDED, NEW SCHMIDT ROUTINE
C 19 NOV 00 - MWS - FIX DIMENSION IN SOSYM COMMON
C 26 OCT 00 - MWS - INTRODUCE MXAO PARAMETER
C 11 JUN 00 - DGF - PRCIVC: ADJUST MXSPIN FOR SINGLE REFERENCE CI
C 19 MAY 00 - DGF - NEW SYMSQT ROUTINE
C  1 MAY 00 - MWS - PRCIVC: CORRECT A FORMAT STATEMENT FOR MANY ORBS.
C 10 JAN 00 - DGF - ADDED KOMVEC ROUTINE
C 21 DEC 99 - DGF - CLEBSCH,IEULER: ADD CLEBSCH-GORDAN COEFS. ROUTINES
C                   INIFAC,ACOS1,DET3,DACOPY,SYMTRZE,ASYMTRZE ADDED.
C                   IRRSTATE,DETFF,SEREP,TRACEW,IDETDEISCO,LRIRREP,
C                   CSFTODET ADDED TO PROVIDE CI EIGENFUNCTION LABELING.
C 12 NOV 98 - GDF - DELETE ISOIN/ISOOUT
C 26 OCT 98 - MAF - TFTRIB: ALLOW FOR USE OF SPHERICAL HARMONICS
C 12 SEP 98 - DGF - PRCIVC: ALL NODES READ DRT FILE, LOOP 980 RANGE
C  6 JAN 98 - DGF - ADD VICLR ROUTINE
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C  2 SEP 97 - MWS - REVERSE TO 80/133 COLUMN OUTPUT FOR NPRINT=7/6
C 16 JUL 97 - GNM - PRATM: CHANGES FOR FRAGMENT ONLY PRINTING
C 14 FEB 97 - MW,SK - PRCIVC: USE VALID INDICES FOR CORE ORBITALS
C  8 JAN 97 - MWS - PRCIVC: CHANGE CORE ORBITAL SKIPPING
C 18 DEC 96 - TLW - ADD MRTRBT ROUTINE
C  1 DEC 96 - MWS - PRCIVC: CORRECT CI CORE OMISSION DURING GRADIENTS
C 26 NOV 96 - VAG - INCREASE CI AO'S TO 768
C 13 JUN 96 - VAG - ADD VARIABLE FOR CI TYPE TO SCFOPT COMMON
C 24 MAY 96 - WC  - PRATM: DECLARE FRGNME AND PTNAM CORRECTLY
C  5 MAR 96 - MWS - ISOIN,ISOOUT: CHANGE SHELL SYMMETRY PACKING
C  9 JAN 96 - MWS - EXTEND FRGMSS COMMON
C 10 OCT 95 - MWS - INCLUDE NEW MRARTR ROUTINE
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 15 JUL 94 - BMB - GAMGEN: ALLOWS DIFFERENT RANGE FOR GRADIENTS
C  3 JUN 94 - PND - PRATM: PRINT FRAGMENT COORDINATES IN DAT FILE
C  5 APR 94 - MWS - INCLUDE MRTRBR; MRARBR,MRTRBR USE DGEMM IF POSSIBLE
C  9 DEC 93 - MWS - TFDS,TRANG,TRANGB: CHANGE DAF RECORD NUMBERS
C                   TFHF,TFHB: PASS DAF RECORD AS ARGUMENT
C 10 AUG 93 - MWS - PUT TFDIAG BACK IN THIS FILE
C  2 AUG 93 - MWS - QMATRX: USE PREVNL ROUTINE FOR DEBUG
C  2 APR 93 - TLW - ADD TRANGB ROUTINE
C  1 APR 93 - PND - PRATM: CHANGES FOR EFFECTIVE FRAGMENTS
C 29 MAR 93 - TLW - TFDS,TFHB,TFHF,TFSQU,TRANG: RECTANGULAR B AND B INV
C 16 JUN 92 - MWS - PARALLELIZE MRARBR,TFSQB,TFTRI; ADD PARALLEL TFSQU;
C                   DELETE MRTRBR,MTATBS
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C 28 FEB 92 - MWS - CHANGE PREVS, DELETE CANVEC,QMAT5,TFDIAG, ADD TFTRIB
C 11 JAN 92 - TLW - PRCIVC: MAKE READS PARALLEL
C 10 JAN 92 - TLW - CHANGE REWIND TO CALL SEQREW
C  7 JAN 92 - TLW - MAKE WRITES PARALLEL; ADD COMMON PAR
C 23 OCT 91 - JHJ - INTR: ADDED RADIUS-CALCULATION FOR ZRF.
C 11 OCT 91 - MWS - DELETED ROUTINE MTARBP, ADDED ROUTINE MRTRBR
C  7 JUL 91 - MWS - ORTHO: Q,S MUST BE PASSED IN, NOT READ FROM DAF.
C                   QMAT5,QMATRX: DO NOT STORE Q MATRIX ON DAF.
C  9 APR 91 - MR  - INTR: NARROW OUTPUT OPTION, MARK SHORT DISTANCES
C 12 SEP 90 - MWS - INTRODUCE MXATM
C  7 AUG 90 - TLW - ADD CF AND CG TO COMMON NSHEL
C  3 JUN 90 - JAM - CHANGE VAX BYTE PACKING TO FULL 8 BITS
C 31 MAY 90 - MWS - ADD ROUTINE PRATM TO PRINT ATOMIC COORDINATES
C 15 DEC 89 - MWS - ADD NEW ROUTINE TFDIAG
C 27 NOV 89 - MWS - BRING ROUTINE TRPOSE TO THIS SOURCE FILE
C 22 SEP 89 - MWS - CHANGE ANGSTROMS-BOHR CONVERSION FACTOR
C 21 SEP 89 - STE - MOVE LOGAND TO MODULE UNPORT
C  6 SEP 89 - MWS - COPY FROM LOWER TRIANGLE OF B IN CPYSQT
C 14 MAR 89 - JAB - MODIFY TRANG TO HANDLE MORE THAN 3N-6 INTERNALS
C 26 FEB 89 - STE - QMAT5,QMATRX: CHANGE TOL=5.0D-08,1.0D-04 TO 1.0D-10
C 27 JAN 89 - MWS - ADD ROUTINE PREVS
C  7 OCT 88 - MWS - INCREASE GUGA AO NUMBER FROM 128 TO 256
C  5 SEP 88 - MWS - ADD NCORBS TO /ORBSET/ AND SKIP MCC'S IN PRCIVC
C 10 AUG 88 - MWS - MXSH,MXGSH,MXGTOT FROM 120,10,440 TO 1000,30,5000
C                   ADD VADD,VCLR,VSUB, RETIRING MODULE VECLIB
C 27 JUL 88 - JAB - ADD EXTRA ARGUMENTS TO TFHF, TFHB FOR VIB. ANALYSIS.
C 19 JUL 88 - MWS - MODIFY PRCIVC TO PRINT THE SPACE ORBITAL PRODUCTS
C  1 JUL 88 - JAM - ADD ROUTINE EXPND
C 30 MAY 88 - MWS - USE PARAMETERS TO DIMENSION COMMON
C  8 MAY 88 - MWS - UP AO-S FROM 255 TO 2047
C 15 APR 88 - MWS - ADD CPYSQT,CPYSS,CPYTSQ,TRPOSQ,MTARBP
C  7 APR 88 - MWS - ADD SLVLEQ, USE DDOT IN TRACEP, CLEAN UP TFSQB
C 29 JAN 88 - MWS - ABORT JOB IF ATOMS ARE TOO CLOSE TOGETHER
C  6 AUG 87 - MWS - INCLUDE ETA VERSION
C  4 MAY 87 - STE - CANVEC: REDUCE WORK SPACE FOR GLDIAG
C 24 APR 87 - STE - TFTRI(TFHB,TFHF),ORTHO: REMOVE IA FROM CALL LIST
C 10 DEC 86 - MWS - PUSQL:  USE COMMA AFTER 1P IN FORMAT (9008)
C  4 NOV 86 - STE - USE PARAMETER STATEMENTS, REMOVE NAV FROM CALLS
C                   PUSQL: CHANGE 9008 TO IMPROVE ORBITAL PRECISION
C 28 OCT 86 - STE - CLARIFY DEPENDENCIES IN MTARBR,MTATBS,SORTIN,TFTRI
C 13 OCT 86 - JAB - ADD NAV ARGS TO DAREAD CALLS IN ORTHO
C 12 AUG 86 - MWS - ADD ROUTINES DGESL, MRARBR, MTARBR, MTATBS
C 22 JUL 86 - MWS - REMOVE RZEROF, USE VCLR INSTEAD
C  9 JUL 86 - JAB - ADD BIT PACKING/UNPACKING FOR CRAY, CELERITY
C  7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
C 27 JUL 86 - MWS - INITIALIZE IN TFHF, COMMENT OUT PREV CONVERSIONS
C 20 JUN 86 - MWS - CLEANUP IN TFHF AND TFHB
C 17 JUN 86 - MWS - LEAVE EIGENVALUES IN PREV IN A.U.
C 10 JAN 86 - LAM - CHANGE TOL IN QMATRX TO 0.5E-7
C  7 NOV 85 - LAM - ACCOMMODATE NMOS.GE.100 IN PUSQL,
C                   PRINT E IN EV IN PREV
C  2 NOV 85 - STE - ONVMGS: FIX LAST VECTOR 0 BUG; USE DDOT,DCOPY
C                   USE GENERIC ABS,EXP,MIN,SQRT; DELETE /MACHIN/
C 10 JUL 85 - MWS - CHANGE IBM INDEX PACKING
C 10 MAR 85 - MWS - INCLUDE LOGAND
C 28 APR 84 - STE - INCLUDE CANVEC,QMAT5 FROM TRFINT
C  9 APR 84 - STE - INCLUDE LINPACK ROUTINES DGEDI,DGEFA
C 23 FEB 84 - STE - INCLUDE PRCIVC
C 17 FEB 84 - STE - INCLUDE ONVMGS,RUNITV,RZEROF,SORTIN
C  1 FEB 84 - STE - INCLUDE PREVNL (FORMERLY XPANVC IN GUGDG)
C 27 DEC 83 - STE - INCLUDE ORDERV, DELETE SCPRD, KILL PRTRI/RUNLAB/
C 29 NOV 83 - STE - FIX FORMAT IN PRTRIL
C 17 NOV 83 - STE - FPS PACKING/UNPACKING IN ISOIN,ISOOUT
C  8 MAY 83 - MWS - PIECES OF OPTX,GAMESS ADDED
C 18 APR 83 - MWS - FIX ERROR IN PRTRIL
C 24 NOV 82 - MWS - PRINT LOWER TRIANGLAR FORM ONLY
C                   SET SMALL VALUES TO ZERO IN TFTRI,TFSQB
C 15 NOV 82 - MWS - FORMAT STATEMENTS IN PRINTOUT ROUTINES
C  4 NOV 82 - MWS - CLEAN ZERO-S IN VECTORS IN QMATRX
C 24 OCT 82 - MWS - INCREASE NUMBER OF COLUMNS IN PRINT ROUTINES
C  1 OCT 82 - MWS - CONVERT TO IBM
C
C*MODULE MTHLIB  *DECK CPYSQT
      SUBROUTINE CPYSQT(B,A,NA,INCA)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(*),B(NA,NA)
      PARAMETER (SMALL=1.0D-10, ZERO=0.0D+00)
C
C     ---- COPY SYMMETRIC, SQUARE B (N X N) TO TRIANGULAR A ----
C     THE INCREMENT BETWEEN ELEMENTS OF A WILL USUALLY BE 1.
C     TINY NUMBERS GET ZEROED OUT.
C
      IJ=1
      DO 200 I=1,NA
         DO 100 J=1,I
            VAL = B(I,J)
            IF(ABS(VAL).LT.SMALL) VAL=ZERO
            A(IJ) = VAL
            IJ = IJ + INCA
  100    CONTINUE
  200 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK CPYTSQ
      SUBROUTINE CPYTSQ(A,B,NA,INCA)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(*),B(NA,NA)
C
C     ---- COPY TRIANGULAR A TO SQUARE B (NA BY NA) ----
C     THE INCREMENT BETWEEN ELEMENTS OF A WILL USUALLY BE 1.
C
      IJ=1
      DO 200 I=1,NA
         DO 100 J=1,I
            B(I,J) = A(IJ)
            B(J,I) = A(IJ)
            IJ = IJ + INCA
  100    CONTINUE
  200 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK DGECO
      SUBROUTINE DGECO(A,LDA,N,IPVT,RCOND,Z)
      INTEGER LDA,N,IPVT(1)
      DOUBLE PRECISION A(LDA,1),Z(1)
      DOUBLE PRECISION RCOND
      DOUBLE PRECISION DDOT,EK,T,WK,WKM
      DOUBLE PRECISION ANORM,S,DASUM,SM,YNORM
      INTEGER INFO,J,K,KB,KP1,L
C
C        SEE COMMENTS IN DGEFA
C
      ANORM = 0.0D+00
      DO 10 J = 1, N
         ANORM = MAX(ANORM,DASUM(N,A(1,J),1))
   10 CONTINUE
      CALL DGEFA(A,LDA,N,IPVT,INFO)
      EK = 1.0D+00
      DO 20 J = 1, N
         Z(J) = 0.0D+00
   20 CONTINUE
      DO 100 K = 1, N
         IF (Z(K) .NE. 0.0D+00) EK = SIGN(EK,-Z(K))
         IF (ABS(EK-Z(K)) .LE. ABS(A(K,K))) GO TO 30
            S = ABS(A(K,K))/ABS(EK-Z(K))
            CALL DSCAL(N,S,Z,1)
            EK = S*EK
   30    CONTINUE
         WK = EK - Z(K)
         WKM = -EK - Z(K)
         S = ABS(WK)
         SM = ABS(WKM)
         IF (A(K,K) .EQ. 0.0D+00) GO TO 40
            WK = WK/A(K,K)
            WKM = WKM/A(K,K)
         GO TO 50
   40    CONTINUE
            WK = 1.0D+00
            WKM = 1.0D+00
   50    CONTINUE
         KP1 = K + 1
         IF (KP1 .GT. N) GO TO 90
            DO 60 J = KP1, N
               SM = SM + ABS(Z(J)+WKM*A(K,J))
               Z(J) = Z(J) + WK*A(K,J)
               S = S + ABS(Z(J))
   60       CONTINUE
            IF (S .GE. SM) GO TO 80
               T = WKM - WK
               WK = WKM
               DO 70 J = KP1, N
                  Z(J) = Z(J) + T*A(K,J)
   70          CONTINUE
   80       CONTINUE
   90    CONTINUE
         Z(K) = WK
  100 CONTINUE
      S = 1.0D+00/DASUM(N,Z,1)
      CALL DSCAL(N,S,Z,1)
C
C     SOLVE TRANS(L)*Y = W
C
      DO 120 KB = 1, N
         K = N + 1 - KB
         IF (K .LT. N) Z(K) = Z(K) + DDOT(N-K,A(K+1,K),1,Z(K+1),1)
         IF (ABS(Z(K)) .LE. 1.0D+00) GO TO 110
            S = 1.0D+00/ABS(Z(K))
            CALL DSCAL(N,S,Z,1)
  110    CONTINUE
         L = IPVT(K)
         T = Z(L)
         Z(L) = Z(K)
         Z(K) = T
  120 CONTINUE
      S = 1.0D+00/DASUM(N,Z,1)
      CALL DSCAL(N,S,Z,1)
C
      YNORM = 1.0D+00
C
C     SOLVE L*V = Y
C
      DO 140 K = 1, N
         L = IPVT(K)
         T = Z(L)
         Z(L) = Z(K)
         Z(K) = T
         IF (K .LT. N) CALL DAXPY(N-K,T,A(K+1,K),1,Z(K+1),1)
         IF (ABS(Z(K)) .LE. 1.0D+00) GO TO 130
            S = 1.0D+00/ABS(Z(K))
            CALL DSCAL(N,S,Z,1)
            YNORM = S*YNORM
  130    CONTINUE
  140 CONTINUE
      S = 1.0D+00/DASUM(N,Z,1)
      CALL DSCAL(N,S,Z,1)
      YNORM = S*YNORM
C
C     SOLVE  U*Z = V
C
      DO 160 KB = 1, N
         K = N + 1 - KB
         IF (ABS(Z(K)) .LE. ABS(A(K,K))) GO TO 150
            S = ABS(A(K,K))/ABS(Z(K))
            CALL DSCAL(N,S,Z,1)
            YNORM = S*YNORM
  150    CONTINUE
         IF (A(K,K) .NE. 0.0D+00) Z(K) = Z(K)/A(K,K)
         IF (A(K,K) .EQ. 0.0D+00) Z(K) = 1.0D+00
         T = -Z(K)
         CALL DAXPY(K-1,T,A(1,K),1,Z(1),1)
  160 CONTINUE
C     MAKE ZNORM = 1.0
      S = 1.0D+00/DASUM(N,Z,1)
      CALL DSCAL(N,S,Z,1)
      YNORM = S*YNORM
C
      IF (ANORM .NE. 0.0D+00) RCOND = YNORM/ANORM
      IF (ANORM .EQ. 0.0D+00) RCOND = 0.0D+00
      RETURN
      END
C*MODULE MTHLIB  *DECK DGEDI
      SUBROUTINE DGEDI(A,LDA,N,IPVT,DET,WORK,JOB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),DET(2),WORK(*),IPVT(*)
C
C     DGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX
C     USING THE FACTORS COMPUTED BY DGECO OR DGEFA.
C
C     ON ENTRY
C
C        A       DOUBLE PRECISION(LDA, N)
C                THE OUTPUT FROM DGECO OR DGEFA.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  A .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C        IPVT    INTEGER(N)
C                THE PIVOT VECTOR FROM DGECO OR DGEFA.
C
C        WORK    DOUBLE PRECISION(N)
C                WORK VECTOR.  CONTENTS DESTROYED.
C
C        JOB     INTEGER
C                = 11   BOTH DETERMINANT AND INVERSE.
C                = 01   INVERSE ONLY.
C                = 10   DETERMINANT ONLY.
C
C     ON RETURN
C
C        A       INVERSE OF ORIGINAL MATRIX IF REQUESTED.
C                OTHERWISE UNCHANGED.
C
C        DET     DOUBLE PRECISION(2)
C                DETERMINANT OF ORIGINAL MATRIX IF REQUESTED.
C                OTHERWISE NOT REFERENCED.
C                DETERMINANT = DET(1) * 10.0**DET(2)
C                WITH  1.0 .LE. ABS(DET(1)) .LT. 10.0
C                OR  DET(1) .EQ. 0.0 .
C
C     ERROR CONDITION
C
C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS
C        A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED.
C        IT WILL NOT OCCUR IF THE ROUTINES ARE CALLED CORRECTLY
C        AND IF DGECO HAS SET RCOND .GT. 0.0 OR DGEFA HAS SET
C        INFO .EQ. 0 .
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     ROUTINES AND FUNCTIONS
C
C     BLAS DAXPY,DSCAL,DSWAP
C     FORTRAN ABS,MOD
C
C     COMPUTE DETERMINANT
C
      IF (JOB/10 .EQ. 0) GO TO 70
         DET(1) = 1.0D+00
         DET(2) = 0.0D+00
         TEN = 10.0D+00
         DO 50 I = 1, N
            IF (IPVT(I) .NE. I) DET(1) = -DET(1)
            DET(1) = A(I,I)*DET(1)
C        ...EXIT
            IF (DET(1) .EQ. 0.0D+00) GO TO 60
   10       IF (ABS(DET(1)) .GE. 1.0D+00) GO TO 20
               DET(1) = TEN*DET(1)
               DET(2) = DET(2) - 1.0D+00
            GO TO 10
   20       CONTINUE
   30       IF (ABS(DET(1)) .LT. TEN) GO TO 40
               DET(1) = DET(1)/TEN
               DET(2) = DET(2) + 1.0D+00
            GO TO 30
   40       CONTINUE
   50    CONTINUE
   60    CONTINUE
   70 CONTINUE
C
C     COMPUTE INVERSE(U)
C
      IF (MOD(JOB,10) .EQ. 0) GO TO 150
         DO 100 K = 1, N
            A(K,K) = 1.0D+00/A(K,K)
            T = -A(K,K)
            CALL DSCAL(K-1,T,A(1,K),1)
            KP1 = K + 1
            IF (N .LT. KP1) GO TO 90
            DO 80 J = KP1, N
               T = A(K,J)
               A(K,J) = 0.0D+00
               CALL DAXPY(K,T,A(1,K),1,A(1,J),1)
   80       CONTINUE
   90       CONTINUE
  100    CONTINUE
C
C        FORM INVERSE(U)*INVERSE(L)
C
         NM1 = N - 1
         IF (NM1 .LT. 1) GO TO 140
         DO 130 KB = 1, NM1
            K = N - KB
            KP1 = K + 1
            DO 110 I = KP1, N
               WORK(I) = A(I,K)
               A(I,K) = 0.0D+00
  110       CONTINUE
            DO 120 J = KP1, N
               T = WORK(J)
               CALL DAXPY(N,T,A(1,J),1,A(1,K),1)
  120       CONTINUE
            L = IPVT(K)
            IF (L .NE. K) CALL DSWAP(N,A(1,K),1,A(1,L),1)
  130    CONTINUE
  140    CONTINUE
  150 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK DGEFA
      SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,*),IPVT(*)
C
C     DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION.
C
C     DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED
C     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
C     (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) .
C
C     ON ENTRY
C
C        A       DOUBLE PRECISION(LDA, N)
C                THE MATRIX TO BE FACTORED.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  A .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C     ON RETURN
C
C        A       AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
C                WHICH WERE USED TO OBTAIN IT.
C                THE FACTORIZATION CAN BE WRITTEN  A = L*U  WHERE
C                L  IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C                TRIANGULAR MATRICES AND  U  IS UPPER TRIANGULAR.
C
C        IPVT    INTEGER(N)
C                AN INTEGER VECTOR OF PIVOT INDICES.
C
C        INFO    INTEGER
C                = 0  NORMAL VALUE.
C                = K  IF  U(K,K) .EQ. 0.0 .  THIS IS NOT AN ERROR
C                     CONDITION FOR THIS ROUTINE, BUT IT DOES
C                     INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO
C                     IF CALLED.  USE  RCOND  IN DGECO FOR A RELIABLE
C                     INDICATION OF SINGULARITY.
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     ROUTINES AND FUNCTIONS
C
C     BLAS DAXPY,DSCAL,IDAMAX
C
C     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
      INFO = 0
      NM1 = N - 1
      IF (NM1 .LT. 1) GO TO 70
      DO 60 K = 1, NM1
         KP1 = K + 1
C
C        FIND L = PIVOT INDEX
C
         L = IDAMAX(N-K+1,A(K,K),1) + K - 1
         IPVT(K) = L
C
C        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
         IF (A(L,K) .EQ. 0.0D+00) GO TO 40
C
C           INTERCHANGE IF NECESSARY
C
            IF (L .EQ. K) GO TO 10
               T = A(L,K)
               A(L,K) = A(K,K)
               A(K,K) = T
   10       CONTINUE
C
C           COMPUTE MULTIPLIERS
C
            T = -1.0D+00/A(K,K)
            CALL DSCAL(N-K,T,A(K+1,K),1)
C
C           ROW ELIMINATION WITH COLUMN INDEXING
C
            DO 30 J = KP1, N
               T = A(L,J)
               IF (L .EQ. K) GO TO 20
                  A(L,J) = A(K,J)
                  A(K,J) = T
   20          CONTINUE
               CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
   30       CONTINUE
         GO TO 50
   40    CONTINUE
            INFO = K
   50    CONTINUE
   60 CONTINUE
   70 CONTINUE
      IPVT(N) = N
      IF (A(N,N) .EQ. 0.0D+00) INFO = N
      RETURN
      END
C*MODULE MTHLIB  *DECK DGESL
      SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,N),B(N),IPVT(N)
C
C     DGESL SOLVES THE DOUBLE PRECISION SYSTEM
C     A * X = B  OR  TRANS(A) * X = B
C     USING THE FACTORS COMPUTED BY DGECO OR DGEFA.
C
C     ON ENTRY
C
C        A       DOUBLE PRECISION(LDA, N)
C                THE OUTPUT FROM DGECO OR DGEFA.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  A .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C        IPVT    INTEGER(N)
C                THE PIVOT VECTOR FROM DGECO OR DGEFA.
C
C        B       DOUBLE PRECISION(N)
C                THE RIGHT HAND SIDE VECTOR.
C
C        JOB     INTEGER
C                = 0         TO SOLVE  A*X = B ,
C                = NONZERO   TO SOLVE  TRANS(A)*X = B  WHERE
C                            TRANS(A)  IS THE TRANSPOSE.
C
C     ON RETURN
C
C        B       THE SOLUTION VECTOR  X .
C
C     ERROR CONDITION
C
C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
C        ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES SINGULARITY
C        BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
C        SETTING OF LDA .  IT WILL NOT OCCUR IF THE ROUTINES ARE
C        CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0
C        OR DGEFA HAS SET INFO .EQ. 0 .
C
C     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX
C     WITH  P  COLUMNS
C           CALL DGECO(A,LDA,N,IPVT,RCOND,Z)
C           IF (RCOND IS TOO SMALL) GO TO ...
C           DO 10 J = 1, P
C              CALL DGESL(A,LDA,N,IPVT,C(1,J),0)
C        10 CONTINUE
C
C     LINPACK. THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     ROUTINES AND FUNCTIONS
C
C     BLAS DAXPY,DDOT
C
      NM1 = N - 1
      IF (JOB .NE. 0) GO TO 50
C
C        JOB = 0 , SOLVE  A * X = B
C        FIRST SOLVE  L*Y = B
C
         IF (NM1 .LT. 1) GO TO 30
         DO 20 K = 1, NM1
            L = IPVT(K)
            T = B(L)
            IF (L .EQ. K) GO TO 10
               B(L) = B(K)
               B(K) = T
   10       CONTINUE
            CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1)
   20    CONTINUE
   30    CONTINUE
C
C        NOW SOLVE  U*X = Y
C
         DO 40 KB = 1, N
            K = N + 1 - KB
            B(K) = B(K)/A(K,K)
            T = -B(K)
            CALL DAXPY(K-1,T,A(1,K),1,B(1),1)
   40    CONTINUE
      GO TO 100
   50 CONTINUE
C
C        JOB = NONZERO, SOLVE  TRANS(A) * X = B
C        FIRST SOLVE  TRANS(U)*Y = B
C
         DO 60 K = 1, N
            T = DDOT(K-1,A(1,K),1,B(1),1)
            B(K) = (B(K) - T)/A(K,K)
   60    CONTINUE
C
C        NOW SOLVE TRANS(L)*X = Y
C
         IF (NM1 .LT. 1) GO TO 90
         DO 80 KB = 1, NM1
            K = N - KB
            B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1)
            L = IPVT(K)
            IF (L .EQ. K) GO TO 70
               T = B(L)
               B(L) = B(K)
               B(K) = T
   70       CONTINUE
   80    CONTINUE
   90    CONTINUE
  100 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK EXPND
      SUBROUTINE EXPND(A,B,N,ISYM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(*),B(N,*)
C
      PARAMETER (ZERO=0.0D+00)
C
C     EXPAND TRIANGULAR MATRIX A TO SQUARE MATRIX B
C     B IS SYMMETRIC IF ISYM=0, ANTISYMMETRIC OTHERWISE
C
      IJ = 0
      IF (ISYM .NE. 0) THEN
         DO 200 J = 1,N
            DO 100 I = 1,J
               IJ = IJ + 1
               B(I,J) = A(IJ)
               B(J,I) = -A(IJ)
  100       CONTINUE
            B(J,J) = ZERO
  200    CONTINUE
      ELSE
         DO 300 J = 1,N
         DO 300 I = 1,J
            IJ = IJ + 1
            B(I,J) = A(IJ)
            B(J,I) = A(IJ)
  300    CONTINUE
      END IF
      RETURN
      END
C*MODULE MTHLIB  *DECK INTR
      SUBROUTINE INTR
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXAO=8192)
      PARAMETER (UNIT=0.52917724924D+00, CLOSE=0.1D+00)
      PARAMETER (RMARK=3.0D+00, ONE=1.00D+00, TWO=2.00D+00)
C
      LOGICAL CROWD,GOPARR,DSKWRK,MASWRK
      INTEGER MARK,STAR,BLANK
C
      DIMENSION R(MXATM)
      DIMENSION MARK(MXATM)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUN,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
C
      DATA STAR /4H *  / , BLANK /4H    /
C
C     ----- CALCULATE ATOM-ATOM DISTANCES -----
C
      IF (MASWRK) WRITE (IW,9008)
      FAC = UNIT
      CROWD = .FALSE.
      INCR = 4
      IF(NPRINT.EQ.6) INCR=7
      MAX = 0
  140 MIN = MAX+1
      MAX = MAX+INCR
      IF (MAX .GT. NAT) MAX = NAT
      IF (MASWRK) THEN
         WRITE (IW,9048)
         WRITE (IW,9068) (ANAM(J),BNAM(J),J = MIN,MAX)
         WRITE (IW,9048)
      END IF
      DO 180 I = 1,NAT
         DO 160 J = MIN,MAX
            RR = (C(1,I)-C(1,J))**2
     *         + (C(2,I)-C(2,J))**2
     *         + (C(3,I)-C(3,J))**2
            RR = SQRT(RR) * FAC
            IF (RR.LE.CLOSE  .AND.  I.NE.J) CROWD=.TRUE.
            IF (RR.GT.RMARK .OR. I.EQ.J) THEN
              MARK(J) = BLANK
            ELSE
              MARK(J) = STAR
            END IF
            R(J) = RR
  160    CONTINUE
         IF (MASWRK) WRITE (IW,9088)
     *      I,ANAM(I),BNAM(I),(R(J),MARK(J),J = MIN,MAX)
  180 CONTINUE
      IF (MAX .LT. NAT) GO TO 140
C
      IF (MASWRK) WRITE(IW,9090)RMARK
      IF(CROWD) THEN
         IF (MASWRK) WRITE(IW,9100) CLOSE
C
C           FMO RUNS ARE ALLOWED TO HAVE ATOMS AT THE SAME PLACE
C           (CROWDED) (DUE TO GHOST ATOMS IN ESPS).
C
         IF(NFG.EQ.0) THEN
           CALL ABRT
           STOP
         ENDIF
      END IF
C
      IF(IZRF*ICALC.EQ.1)THEN
        CALL ACALC(RADIUS)
        IF (MASWRK) WRITE(IW,9110)RADIUS
        RADIUS = RADIUS/UNIT
        GZRF=TWO*(DIELEC-ONE)/((TWO*DIELEC+ONE)*RADIUS**3)
        EBORN=(1-DIELEC)*ICH**2/(TWO*DIELEC*RADIUS)
      END IF
C
      RETURN
C
 9008 FORMAT(/,10X,30HINTERNUCLEAR DISTANCES (ANGS.),/,10X,30(1H-))
 9048 FORMAT(1X)
 9068 FORMAT(17X,7(3X,A8,A2,2X))
 9088 FORMAT(I3,2X,A8,A2,F15.7,A4,6(F11.7,A4))
 9090 FORMAT(/'  * ... LESS THAN',F7.3/)
 9100 FORMAT(1X,'**** THERE ARE ATOMS LESS THAN',F8.3,
     *          ' APART, QUITTING... ****')
 9110 FORMAT(/,1X,'ADJUSTED CAVITY RADIUS = ',F5.2)
      END
C*MODULE MTHLIB  *DECK LOCAL
      SUBROUTINE LOCAL(X,Y,Z,XS,YS,ZS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      COMMON /FRAME / U1,U2,U3,V1,V2,V3,W1,W2,W3,X0,Y0,Z0
C
C     CALCULATE THE COORDINATES (XS,YS,ZS) OF A POINT IN THE LOCAL
C     FRAME GIVEN THE COORDINATES (X,Y,Z) IN THE MASTER FRAME
C
      XS = U1*(X-X0)+U2*(Y-Y0)+U3*(Z-Z0)
      YS = V1*(X-X0)+V2*(Y-Y0)+V3*(Z-Z0)
      ZS = W1*(X-X0)+W2*(Y-Y0)+W3*(Z-Z0)
      RETURN
      END
C*MODULE MTHLIB  *DECK LSTRNG
      INTEGER FUNCTION LSTRNG(STRING,LEN)
      CHARACTER*(*) STRING
      DO I=1,LEN
         IF(STRING(I:I).EQ.' ') THEN
            LSTRNG=I-1
            RETURN
         END IF
      ENDDO
      LSTRNG=LEN
      RETURN
      END
C*MODULE MTHLIB  *DECK MRARBR
      SUBROUTINE MRARBR(A,LDA,NA,MA,B,LDB,MB,C,LDC)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION A(LDA,MA),B(LDB,MB),C(LDC,MB)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00)
*BL3  PARAMETER (ONE=1.0D+00)
C
C     ----- SEQUENTIAL MATRIX MULTIPLY: C = A * B -----
C        A    - THE INPUT REAL NA X MA RECTANGULAR MATRIX
C        B    - THE INPUT REAL MA X MB RECTANGULAR MATRIX
C        C    - THE OUTPUT PRODUCT NA X MB MATRIX
C        LDA  - THE LEADING DIMENSION OF MATRIX A
C        NA   - THE EFFECTIVE ROW DIMENSION OF MATRICES A AND C
C        MA   - THE COLUMN DIMENSION OF MATRIX A,
C               AND EFFECTIVE ROW DIMENSION OF MATRIX B
C        LDB  - THE LEADING DIMENSION OF MATRIX B
C        MB   - THE COLUMN DIMENSION OF MATRICES B AND C
C        LDC  - THE LEADING DIMENSION OF MATRIX C
C     AUTHOR = STEVE ELBERT, 31 OCT 1979
C
      IF(LDA.LT.NA .OR. LDB.LT.MA .OR. LDC.LT.NA) GO TO 800
C
*BL3  CALL DGEMM('N','N',NA,MB,MA,ONE,A,LDA,B,LDB,ZERO,C,LDC)
*BL3  IF(ONE.NE.ZERO) RETURN
C
      M=MA
      IF(MOD(M,2).NE.0) M=M-1
C
      DO 300 I=1,NA
         DO 200 J=1,MB
            CIJ=ZERO
            IF(M.NE.MA) CIJ=A(I,MA)*B(MA,J)
            IF(MA.GT.1) THEN
               DO 100 K=1,M,2
                  CIJ=CIJ + A(I,K)*B(K,J) + A(I,K+1)*B(K+1,J)
  100          CONTINUE
            END IF
            C(I,J)=CIJ
  200    CONTINUE
  300 CONTINUE
      RETURN
C
  800 CONTINUE
      IF(MASWRK) WRITE(IW,900) LDA,NA,MA,LDB,MB,LDC
      CALL ABRT
      STOP
C
  900 FORMAT(/1X,'ERROR IN CALL TO -MRARBR-'/
     *        1X,'LDA,NA,MA,LDB,MB,LDC=',6I10)
      END
C*MODULE MTHLIB  *DECK MRARTR
      SUBROUTINE MRARTR(A,LDA,LL,MM,B,LDB,NN,C,LDC)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION A(LDA,MM),B(LDB,MM),C(LDC,NN)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00)
*BL3  PARAMETER (ONE=1.0D+00)
C
C     ----- SEQUENTIAL MATRIX MULTIPLY: C = A * B-TRANSPOSE -----
C        A    - THE INPUT REAL LL X MM RECTANGULAR MATRIX
C        B    - THE INPUT REAL NN X MM RECTANGULAR MATRIX
C        C    - THE OUTPUT PRODUCT LL X NN MATRIX
C        LDA  - THE LEADING DIMENSION OF MATRIX A
C        LL   - THE EFFECTIVE ROW DIMENSION OF MATRICES A AND C
C        MM   - THE COLUMN DIMENSION OF MATRIX A,
C               AND EFFECTIVE ROW DIMENSION OF MATRIX B
C        LDB  - THE LEADING DIMENSION OF MATRIX B
C        NN   - THE COLUMN DIMENSION OF MATRICES B AND C
C        LDC  - THE LEADING DIMENSION OF MATRIX C
C     THE ARGUMENTS MAP DIRECTLY TO A LEVEL 3 BLAS CALL.
C
      IF(LDA.LT.LL .OR. LDB.LT.NN .OR. LDC.LT.LL) GO TO 800
C
*BL3  CALL DGEMM('N','T',LL,NN,MM,ONE,A,LDA,B,LDB,ZERO,C,LDC)
*BL3  IF(ONE.NE.ZERO) RETURN
C
      KMAX=MM
      IF(MOD(KMAX,2).NE.0) KMAX=KMAX-1
C
      DO 300 I=1,LL
         DO 200 J=1,NN
            CIJ=ZERO
            IF(KMAX.NE.MM) CIJ=A(I,MM)*B(J,MM)
            IF(MM.GT.1) THEN
               DO 100 K=1,KMAX,2
                  CIJ=CIJ + A(I,K)*B(J,K) + A(I,K+1)*B(J,K+1)
  100          CONTINUE
            END IF
            C(I,J)=CIJ
  200    CONTINUE
  300 CONTINUE
      RETURN
C
  800 CONTINUE
      IF(MASWRK) WRITE(IW,900) LDA,LL,MM,LDB,NN,LDC
      CALL ABRT
      STOP
C
  900 FORMAT(/1X,'ERROR IN CALL TO -MRARTR-'/
     *        1X,'LDA,LL,MM,LDB,NN,LDC=',6I10)
      END
C*MODULE MTHLIB  *DECK MRTRBR
      SUBROUTINE MRTRBR(A,LDA,NA,MA,B,LDB,MB,C,LDC)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(LDA,MA),B(LDB,MB),C(LDC,MB)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
C
      PARAMETER (ZERO=0.0D+00)
*BL3  PARAMETER (ONE=1.0D+00)
C
C     ----- MATRIX MULTIPLY: C = A-TRANSPOSE * B -----
C        A    - THE INPUT REAL NA BY MA RECTANGULAR MATRIX
C        B    - THE INPUT REAL NA BY MB RECTANGULAR MATRIX
C        C    - THE OUTPUT PRODUCT MA BY MB MATRIX
C        LDA  - THE LEADING DIMENSION OF MATRIX A
C        NA   - THE EFFECTIVE ROW DIMENSION OF MATRICES A AND B
C        MA   - THE COLUMN DIMENSION OF MATRIX A
C               AND EFFECTIVE ROW DIMENSION OF MATRIX C
C        LDB  - THE LEADING DIMENSION OF MATRIX B
C        MB   - THE COLUMN DIMENSION OF MATRICES B AND C
C        LDC  - THE LEADING DIMENSION OF MATRIX C
C     AUTHOR = STEVE ELBERT, 31 OCT 1979
C
      IF(LDA.LT.NA .OR. LDB.LT.NA .OR. LDC.LT.MA) GO TO 800
C
*BL3  CALL DGEMM('T','N',MA,MB,NA,ONE,A,LDA,B,LDB,ZERO,C,LDC)
*BL3  IF(ONE.NE.ZERO) RETURN
C
      DO 110 J=1,MB
         DO 100 I=1,MA
            C(I,J)=ZERO
  100    CONTINUE
  110 CONTINUE
C
      DO 140 J=1,MB
         DO 130 K=1,NA
            BKJ=B(K,J)
            IF(BKJ.EQ.ZERO) GO TO 130
               DO 120 I=1,MA
                  C(I,J)=C(I,J)+A(K,I)*BKJ
  120          CONTINUE
  130    CONTINUE
  140 CONTINUE
      RETURN
C
  800 WRITE(IW,900) LDA,NA,MA,LDB,MB,LDC
      CALL ABRT
      STOP
C
  900 FORMAT(1X,'ERROR CALLING -MRTRBR-'/
     *       1X,'LDA,NA,MA,LDB,MB,LDC=',6I10)
      END
C*MODULE MTHLIB  *DECK MRTRBT
      SUBROUTINE MRTRBT(A,NAA,NA,MA,B,NB,AB,J1,I1)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DOUBLE PRECISION A(NAA,MA),B(NB,MA),AB(1)
C
      PARAMETER (ZERO=0.0D+00)
C
C* 31 OCT 1979
C*
C*FUNCTION    - TO MULTIPLY TRANSPOSE OF RECTANGULAR MATRIX A
C*              TIMES RECTANGULAR MATRIX B AND GET SYMMETRIC MATRIX AB
C*
C*PARAMETERS
C*       A    - THE INPUT REAL NA BY MA RECTANGULAR MATRIX
C*       B    - THE INPUT REAL NA BY MA RECTANGULAR MATRIX
C*       AB   - THE OUTPUT SYMMETRIC MATRIX AB OF ORDER MA
C*       NAA  - THE INPUT ROW DIMENSION OF MATRIX A
C*       NA   - THE EFFECTIVE ROW DIMENSION OF MATRICES A AND B
C*       MA   - THE COLUMN DIMENSION OF MATRICES A AND B
C*              AND ORDER OF MATRIX AB
C*       NB   - THE INPUT ROW DIMENSION OF MATRIX B
C*       I1,J1 - INDICES OF FIRST ELEMENT OF OUTPUT MATRIX AB
C
  900 FORMAT(13H0NAA,NA,MA,NB,6I10)
  901 FORMAT(9H0MATRIX A/(1X,10F13.6))
  902 FORMAT(9H0MATRIX B/(1X,10F13.6))
C
C-----------------------------------------------------------------------
C
      I0=I1
      MM=(MA*MA+MA)/2 - (J1*J1-J1)/2 - I0 + 1
      IF(NAA.LT.NA .OR. NB.LT.NA) GO TO 800
         DO 110 I=1,MM
            AB(I)=ZERO
  110    CONTINUE
         II=0
         DO 140 J=J1,MA
            DO 130 K=1,NA
               BKJ=B(K,J)
               IF(BKJ.EQ.ZERO) GO TO 130
                  IJ=II
                  DO 120 I=I0,J
                     IJ=IJ+1
                     AB(IJ)=AB(IJ)+A(K,I)*BKJ
  120             CONTINUE
  130       CONTINUE
            II=II+J-I0+1
            I0=1
  140    CONTINUE
         RETURN
C
  800 WRITE(6,900) NAA,NA,MA,NB
      WRITE(6,901) A
      WRITE(6,902) B
      CALL ABRT
      STOP
      END
C*MODULE MTHLIB  *DECK MTARBR
      SUBROUTINE MTARBR(A,NA,B,MB,AB,NAB,INCA)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(*),B(NA,MB),AB(NAB,MB)
C
      PARAMETER (ZERO=0.0D+00)
C
C* 31 OCT 1979
C*
C*FUNCTION   - TO MULTIPLY SYMMETRIC MATRIX A
C*             TIMES RECTANGULAR MATRIX B AND GET RECTANGULAR MATRIX AB
C*
C*PARAMETERS
C*    A     - THE INPUT REAL SYMMETRIC MATRIX OF ORDER NA
C*            STORED IN SYMMETRIC STOAGE MODE.
C*    B     - THE INPUT REAL NA BY MB RECTANGULAR MATRIX
C*    NA    - THE ORDER OF MATRIX A
C*    MB    - THE COLUMN DIMENSION OF MATRICES B AND AB
C*    AB    - THE OUTPUT PRODUCT NA BY MB MATRIX
C*    NAB   - THE INPUT ROW DIMENSION OF MATRIX AB
C*    INCA  - ADDRESS DIFFERENCE OF TWO ADJACENT ELEMENTS OF A
C*
      INC=INCA
C
C        PROCESS DIAGONAL ELEMENTS OF INPUT MATRIX A
C
      IJ=1-INC
      DO 120 I=1,NA
         IJ=IJ+I*INC
         AIJ=A(IJ)
         DO 110 K=1,MB
            AB(I,K)=AIJ*B(I,K)
  110    CONTINUE
  120 CONTINUE
      IF(NA.EQ.1) RETURN
C
C        PROCESS OFF-DIAGONAL ELEMENTS OF INPUT MATRIX A
C
      IJ=1-INC
      DO 150 I=2,NA
         IJ=IJ+INC
         IM1=I-1
         DO 140 J=1,IM1
            IJ=IJ+INC
            AIJ=A(IJ)
            IF(AIJ.EQ.ZERO) GO TO 140
               CALL DAXPY(MB,AIJ,B(I,1),NA,AB(J,1),NAB)
               CALL DAXPY(MB,AIJ,B(J,1),NA,AB(I,1),NAB)
  140    CONTINUE
  150 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK ONVMGS
      SUBROUTINE ONVMGS(NROWS,NCOLS,VECT,NDROP)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VECT(NROWS,NCOLS)
      PARAMETER (ONE=1.0D+00)
C
C* 23 MAY 1982 - STE
C*
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY - USDOE)  DATE: MAR 1983
C*       (BASED ON NOTES BY CLEVE MOLER PRESENTED AT THE NRCC WORKSHOP
C*        "NUMERICAL ALGORITHMS IN CHEMISTRY:  ALGEBRAIC METHODS")
C*
C*    PURPOSE -
C*       ORTHONORALIZE VECTORS USING THE MODIFIED GRAM-SCHMIDT METHOD
C*       *    *        *                 *        *    *
C*
C*    ON ENTRY -
C*       NROWS  - INTEGER
C*                ROW DIMENSION OF VECTOORS
C*       NCOLS  - INTEGER
C*                COLUMN DIMENSION OF VECTORS
C*       VECT   - W.P. REAL (NROWS,NCOLS)
C*                VECTORS TO BE ORTHONORMALIZED
C*
C*    ON EXIT -
C*       VEXT   - W.P. REAL (NROWS,NCOLS)
C*                ORTHONORMAL SET OF VECTORS
C*       NDROP  - INTEGER
C*                NO. OF LINEARLY DEPENDENT VECTORS
C*                (DEPENDENT VECTORS ARE SET TO ZERO)
C*
C*    COMPLEXITY -
C*       NROWS * NCOLS**2
C*
      NDROP=0
      IF(NCOLS.EQ.1) GO TO 130
         NCOLM1=NCOLS-1
         DO 120 I=1,NCOLM1
            VLENSQ=DDOT(NROWS,VECT(1,I),1,VECT(1,I),1)
            IF(ONE+VLENSQ.NE.ONE) GO TO 100
C
C                 VECTOR IS LINEARLY DEPENDENT
C
               NDROP=NDROP+1
               CALL VCLR(VECT(1,I),1,NROWS)
               GO TO 120
C
  100       CONTINUE
            VNORM=ONE/SQRT(VLENSQ)
            CALL DSCAL(NROWS,VNORM,VECT(1,I),1)
            IPLUS1 = I + 1
            DO 110 J=IPLUS1,NCOLS
               VOVLP=DDOT(NROWS,VECT(1,I),1,VECT(1,J),1)
               CALL DAXPY(NROWS,-VOVLP,VECT(1,I),1,VECT(1,J),1)
  110       CONTINUE
  120    CONTINUE
  130 CONTINUE
C
C        DO LAST VECTOR
C
      VLENSQ=DDOT(NROWS,VECT(1,NCOLS),1,VECT(1,NCOLS),1)
      IF(ONE+VLENSQ.NE.ONE) GO TO 140
C
C           LAST VECTOR IS LINEARLY DEPENDENT
C
         NDROP=NDROP+1
         CALL VCLR(VECT(1,NCOLS),1,NROWS)
         RETURN
C
  140 CONTINUE
      VNORM=ONE/SQRT(VLENSQ)
      CALL DSCAL(NROWS,VNORM,VECT(1,NCOLS),1)
      RETURN
      END
C*MODULE MTHLIB  *DECK ORDERV
      SUBROUTINE ORDERV(C,D,N,X,IA,IB,NBASIS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION C(*),D(IA),N(IA),X(IA)
C
C          VERSION 16 - SEPTEMBER 1972 - S. ELBERT
C        ORDER THE IB LOWEST EIGENVALUES D AND EIGENVECTORS C
C        FROM A SET OF IA EIGENVECTORS OF LENGTH NBASIS
C        N AND X ARE SCRATCH REGIONS OF LENGTH IA AND NBASIS
C
C        INITIALIZE PERMUTATION VECTOR
C
      NA=NBASIS
      DO 110 I=1,IA
         N(I)=I
  110 CONTINUE
      IN=IA-1
      IF(IN.LE.0) RETURN
C
C        PUT ELEMENTS OF D IN ASCENDING ORDER
C
      DO 130 I=1,IN
         DO 120 J=1,I
            K=I-J+1
            IF(D(K+1).GE.D(K))GO TO 130
               X(1)=D(K)
               D(K)=D(K+1)
               D(K+1)=X(1)
               JA=N(K)
               N(K)=N(K+1)
               N(K+1)=JA
  120    CONTINUE
  130 CONTINUE
      IF(IB+1.GE.IA) GO TO 160
C
C        RESTORE ORIGINAL ORDER FOR ELEMENTS IB+1 TO IA
C
      IC=IB+1
      DO 150 I=IC,IN
         DO 140 J=IC,I
            K=I-J+IC
            IF(N(K+1).GE.N(K)) GO TO 150
C
C                    SWAP POSITIONS
C
               X(1)=D(K)
               D(K)=D(K+1)
               D(K+1)=X(1)
               JA=N(K)
               N(K)=N(K+1)
               N(K+1)=JA
  140    CONTINUE
  150 CONTINUE
  160 CONTINUE
C
C          REARRANGE EIGENVECTORS
C
      DO 250 KH=1,IA
         K=N(KH)
         IF(KH.EQ.K) GO TO 250
         IF(K.LE.0) GO TO 250
            IG=NA*K
            IH=IG
            DO 210 I=1,NA
               X(I)=C(IH)
               IH=IH-1
  210       CONTINUE
  220       CONTINUE
               KL=N(K)
               N(K)=-N(K)
               K=KL
               IH=NA*K
               IN=IG
               IG=IH
               IN = IN - NA
               IH = IH - NA
               CALL DCOPY(NA,C(IH+1),1,C(IN+1),1)
            IF(K.NE.KH) GO TO 220
            DO 240 I=1,NA
               C(IG)=X(I)
               IG=IG-1
  240       CONTINUE
  250 CONTINUE
      DO 260 I=1,IA
         N(I) = ABS(N(I))
  260 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK ORTHO
      SUBROUTINE ORTHO(Q,S,V,T,N,L0,L1,L2,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION Q(NDIM,L0),S(L2),V(NDIM,L0),T(L0)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- ORTHOGONALIZE ORBITALS -----
C     OVERLAP -S- AND CANONICAL ORTHONORMAL ORBITALS -Q-
C     MUST BE GIVEN ON ENTRY, AND WILL BE UNCHANGED.
C     THE FIRST -N- VECTORS IN -V- WILL BE GRAM-SCHMIDT
C     ORTHOGONALIZED, THEN L0-N VECTORS WILL BE GENERATED
C     AS THE ORTHOGONAL COMPLEMENT TO THE FIRST -N-.
C     -T- IS WORKING STORAGE ONLY.
C
      DO 180 J = 1,N
         IK = 0
         DO 140 I = 1,L1
            IM1 = I-1
            DUM = ZERO
            VDUM = V(I,J)
            IF(IM1.GT.0) THEN
               DO 100 K = 1,IM1
                  IK = IK+1
                  DUM = DUM+S(IK)*V(K,J)
                  T(K) = T(K)+S(IK)*VDUM
  100          CONTINUE
            END IF
            IK = IK+1
            T(I) = DUM+S(IK)*VDUM
  140    CONTINUE
         DO 170 I = 1,L1
            V(I,J) = T(I)
  170    CONTINUE
  180 CONTINUE
C
      DO 440 J = 1,N
         DO 410 I = 1,L0
            T(I) = DDOT(L1,Q(1,I),1,V(1,J),1)
  410    CONTINUE
         DO 430 I = 1,L0
            V(I,J) = T(I)
  430    CONTINUE
  440 CONTINUE
C
      CALL SCHMD(V,N,L0,NDIM)
C
      RETURN
      END
C*MODULE MTHLIB  *DECK PRATM
      SUBROUTINE PRATM(LU,IMODE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXSH=5000, MXATM=2000, MXAO=8192, MXPT=2000,
     *           MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      LOGICAL SKIP(MXATM)
      LOGICAL FGONLY
C
      CHARACTER*4 UNITS
      CHARACTER*6 FRGNAM
      CHARACTER*8 FRGNME,PTNAM
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
      COMMON /FRGMSS/ FPMASS(MXPT),FMASS(MXFGPT),
     *                FPNUC(MXPT),FGNUC(MXFGPT)
      COMMON /FRAGMT/ XCRD(3,MXFRG),YCRD(3,MXFRG),ZCRD(3,MXFRG),
     *                PTNAM(3,MXFRG),FRGNAM(MXFRG)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
      PARAMETER (TOANGS=0.52917724924D+00)
      PARAMETER (ZERO=0.0D+00)
C
      FGONLY = NAT.EQ.0  .AND.  NUM.EQ.0
C
C     ----- WRITE ATOMIC COORDINATES TO UNIT -LU- -----
C     IF ABS(IMODE).EQ.1, ALL ATOMS
C     IF ABS(IMODE).EQ.2, ONLY THE SYMMETRY UNIQUE ATOMS
C     PRINT IN BOHR IF IMODE IS NEGATIVE, ANGSTROMS OTHERWISE.
C
      MODE = IABS(IMODE)
      IF(MODE.EQ.0) RETURN
C
      IF(IMODE.GT.0) THEN
         UNIT = TOANGS
         UNITS='ANGS'
      ELSE
         UNIT = 1.0D+00
         UNITS='BOHR'
      END IF
C
      IF (FGONLY) THEN
         IF(MASWRK) WRITE (LU,*) ' '
         GO TO 300
      END IF
      IF(MODE.EQ.1 ) THEN
         IF (MASWRK) WRITE(LU,9000) UNITS
         DO 100 IAT = 1,NAT
            SKIP(IAT) = .FALSE.
  100    CONTINUE
         GO TO 300
      ELSE
         IF (MASWRK) WRITE(LU,9010) UNITS
         DO 110 IAT = 1,NAT
            SKIP(IAT) = .TRUE.
  110    CONTINUE
      END IF
C
      DO 260 IAT=1,NAT
         DO 240 IT = 1,NT
            IC = MAPCTR(IAT,IT)
            IF (IC .GT. IAT) GO TO 260
  240    CONTINUE
         SKIP(IAT) = .FALSE.
  260 CONTINUE
C
  300 CONTINUE
      NFRN=0
      IF(NFRG.GT.0)THEN
         KMASS=0
         DO 302 IFRG=1,NFRG
         DO 303 III=1,NMPTS(IFRG)
         IF(FMASS(III+KMASS).GT.ZERO) NFRN=NFRN+1
 303     CONTINUE
         KMASS=KMASS+NMPTS(IFRG)
 302     CONTINUE
      END IF
C
      IF (MASWRK  .AND.  .NOT.FGONLY) WRITE(LU,9020)
C
      DO 310 IAT=1,NAT
         IF(SKIP(IAT)) GO TO 310
         ZNUC = ZAN(IAT) + IZCORE(IAT)
         X = UNIT * C(1,IAT)
         Y = UNIT * C(2,IAT)
         Z = UNIT * C(3,IAT)
         IF (MASWRK) WRITE(LU,9030) ANAM(IAT),BNAM(IAT),ZNUC,X,Y,Z
  310 CONTINUE
      IF(NFRG.GT.0) THEN
      IF (MASWRK) WRITE(LU,9011) UNITS
      IF (MASWRK) WRITE(LU,9021)
C     DO 311 IAT=1,NMTTPT
      IAT=0
      KMASS=0
      DO 311 IFRG=1,NFRG
      IF(MASWRK) WRITE(LU,9022)FRGNAM(IFRG)
      DO 312 III=1,NMPTS(IFRG)
      IAT=IAT+1
      IF(FMASS(III+KMASS).GT.ZERO) THEN
         X = UNIT * EFC(1,IAT)
         Y = UNIT * EFC(2,IAT)
         Z = UNIT * EFC(3,IAT)
         IF (MASWRK) WRITE(LU,9031) FRGNME(IAT),X,Y,Z
      END IF
  312 CONTINUE
      KMASS=KMASS+NMPTS(IFRG)
  311 CONTINUE
C     IAT=0
C     DO 911 IFRG=1,NFRG
C     DO 911 III=1,NMPTS(IFRG)
C        IAT=IAT+1
C      IF(FMASS(III,IFRG).GT.ZERO)THEN
C        X = UNIT * EFC(1,IAT)
C        Y = UNIT * EFC(2,IAT)
C        Z = UNIT * EFC(3,IAT)
C      END IF
C 911 CONTINUE
      END IF
      RETURN
C
 9000 FORMAT(1X,'COORDINATES OF ALL ATOMS ARE (',A4,')')
 9010 FORMAT(1X,'COORDINATES OF SYMMETRY UNIQUE ATOMS (',A4,')')
 9011 FORMAT(1X,'COORDINATES OF FRAGMENT MULTIPOLE CENTERS (',A4,')')
 9020 FORMAT(3X,'ATOM',3X,'CHARGE',7X,'X',14X,'Y',14X,'Z'/
     *       1X,60(1H-))
 9021 FORMAT(1X,'MULTIPOLE NAME',8X,'X',14X,'Y',14X,'Z'/
     *       1X,60(1H-))
 9022 FORMAT(1X,'FRAGNAME=',A6)
 9030 FORMAT(1X,A8,A2,F5.1,3F15.10)
 9031 FORMAT(1X,A8,7X,3F15.10)
      END
C*MODULE MTHLIB  *DECK PRCIVC
      SUBROUTINE PRCIVC(NFT11,NSTATE,NWKS,EIGVAL,EIGVEC,PRTTOL,
     *                  LEVIR,ICASE,IECONF,IARC,NOMX,NR4,NROWS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,RUNTRM,ABELPT,SVDSKW,NONABEL,SKIPS,
     *        NONZ,BAIL
C
      PARAMETER (MXAO=8192, MAXCP=4096, MAXL=3, MXATM=2000,
     *           MXSH=5000, MAXIRR=14)
C
      DIMENSION EIGVEC(NWKS,NSTATE), EIGVAL(NSTATE)
      DIMENSION LEVIR(NOMX),ICASE(NOMX),IECONF(NOMX),IARC(NR4)
      DIMENSION ISHIFT(4),WEIGHT(MAXIRR)
C
      COMMON /FMCOM / X(1)
      COMMON /GUGWFN/ NFZC,NMCC,NDOC,NAOS,NBOS,NALP,NVAL,NEXT,NFZV,
     *                IFORS,IEXCIT,ICICI,NOIRR
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,MA,MB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /ORBSET/ NORBMX,NORBS,NCORBS,NLEVS,NA,NB,NC,NSYM,MSYM,
     *                IDOCC,IVAL,IMCC,ISYM(MXAO),ICODE(MXAO),
     *                NLCS(MXAO),LEVPT(MXAO),LEVNR(MXAO),IOUT(MXAO),
     *                NREFS,IEXCT,NFOCI,INTACT
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMREP/ IRPNAM(14),IPA(14),LAMBDA(14),LAMBD0(14),
     *                IADDR1(14),IADDR2(14),IADDR3(14)
      COMMON /SOGUG/  CP(MAXCP),NUNIQ,LIOBP,ISODAF,NSODA,LSTREC,JSODAF,
     *                JSODA,NRECJ
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,NSAFMO
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TRNMOM/ OPERR,SYMTOL,NOCC,NUMCI,MFZC,NUMVEC,ICI,MXRT,NSTAT
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA TRANST/8HTRANSITN/,GUGA/8HGUGA    /
C
C     ----- PRINT COEFFICIENTS OF CI VECTORS LARGER THAN PRTTOL -----
C     THE ELECTRON OCCUPANCY OF THE PARTIALLY OCCUPIED ORBITALS
C     IS ALSO PRINTED OUT.
C
C     MCSCF WAVEFUNCTIONS, ICICI=0,
C     AND CI GRADIENT RUNS, ICICI=1,
C     SHOULD OMIT THE PRINTING OF THE CORE ORBITALS AT THIS POINT.
C
      RUNTRM=RUNTYP.EQ.TRANST
      NAO = NDOC + NAOS + NBOS + NALP + NVAL + NEXT
      NAO1 = NDOC + NAOS + NBOS + NALP + NVAL
      NAEL= NDOC*2 + NAOS + NBOS + NALP
C     TOL=1.0D-04
      SVDSKW=DSKWRK
      NONABEL=.NOT.ABELPT()
      CALL VALFM(LOADFM)
      IF(NONABEL) THEN
         MAXDET=NAO1*NWKS
         MXSPIN=MIN(NAO1+NFOCI,NAEL)
C        SINGLE REFERENCE CI
         IF(IFORS.EQ.0.AND.NFOCI.EQ.0) MXSPIN=NAOS+NBOS+NALP+IEXCIT*2
         IF(MXSPIN.GE.64/NWDVAR) THEN
            WRITE(IW,*) 'MXSPIN ',MXSPIN,'.GE.',64/NWDVAR
            CALL ABRT
         END IF
         MXPRM=2**MXSPIN
         LIGCSF  = LOADFM  + 1
         LCCI    = LIGCSF  + (NAEL*MAXDET-1)/NWDVAR+1
         LFF     = LCCI    + MAXDET
         LCHARR  = LFF     + MAXL*NT*NAO
         LIRRID  = LCHARR  + NT*NIRRED
         LIRROW  = LIRRID  + (NAO-1)/NWDVAR+1
         LNATOMAO= LIRROW  + (NAO-1)/NWDVAR+1
         LCLCAO  = LNATOMAO+ (NAT-1)/NWDVAR+1
         LS      = LCLCAO  + NUM*NUM
         LWORK1  = LS      + (NUM*NUM+NUM)/2
         LFUNSYM = LWORK1  + NUM
         LPROJ   = LFUNSYM + NUM*MAXL
         LPK     = LPROJ   + NUM*MAXL*MAXL
         LIPRIM  = LPK    + MXSPIN*MXPRM
         LCGC    = LIPRIM + (MXSPIN*MXPRM-1)/NWDVAR+1
         LBK     = LCGC   + MXPRM
         LNTRAP  = LBK    + MXSPIN
         LIICAS  = LNTRAP + (MXSPIN-1)/NWDVAR+1
         LAST    = LIICAS + (MXSPIN-1)/NWDVAR+1
      ELSE
         LAST    = LOADFM  + 1
      END IF
      SKIPS=.TRUE.
      IF(CITYP.EQ.GUGA) SKIPS=NOIRR.GT.0
      IF(NONABEL) THEN
         CALL GOTFM(NGOT)
         NEED=LAST-LOADFM-1
         IF(NEED.GT.NGOT) THEN
            IF(MASWRK) WRITE(IW,*) 'SYMMETRY NEEDS MORE MEMORY',NEED,
     *                             ' GOT',NGOT
            LAST = LOADFM  + 1
            SKIPS=.TRUE.
         END IF
      END IF
      IF(RUNTRM) THEN
         SMALL=SYMTOL
      ELSE
         SMALL=1.0D-04
      END IF
      IF(NOIRR.GE.-1) SMALLCSF=SMALL/NT
      IF(NOIRR.LT.-1) SMALLCSF=10**NOIRR
      IF(RUNTRM) THEN
         LCSFIR = LAST
         LSTSYM = LCSFIR + (NWKS-1)/NWDVAR+1
         LAST   = LSTSYM + (NSTATE-1)/NWDVAR+1
      END IF
      IF(RUNTRM.OR.NONABEL) THEN
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
C        SOME OF THE REQUESTS MAY BE OF 0 LENGTH. IS THAT A PROBLEM?
      END IF
      IF(RUNTRM.AND.ABELPT()) THEN
         DSKWRK=.TRUE.
         CALL RAREAD(JSODAF,X(JSODA),X(LCSFIR),(NWKS-1)/NWDVAR+1,
     *               (ICI-1)*NRECJ+5,1)
         DSKWRK=SVDSKW
      END IF
      DSKWRK=SVDSKW
      CALL DERCHK(NDER)
      NCORE = NCORBS
      NPART = NORBS
      IF(ICICI.EQ.1  .AND.  NDER.EQ.1) NPART = NORBS-NCORBS
      IF(ICICI.EQ.0) THEN
         NCORE = NMCC
         NPART = NORBS - NCORE
      END IF
C
      IF(NONABEL.AND..NOT.SKIPS) THEN
         IF(MASWRK) WRITE(IW,9200) NEED
         CALL DAREAD(IDAF,IODA,X(LCLCAO),NUM*NUM,15,0)
         CALL DAREAD(IDAF,IODA,X(LS),(NUM*NUM+NUM)/2,12,0)
         BAIL = .TRUE.
         CALL SEREP(NAO,NCORE,X(LFF),X(LCHARR),X(LIRRID),X(LIRROW),
     *              X(LCLCAO),X(LS),X(LWORK1),X(LFUNSYM),X(LPROJ),
     *              SMALL,.FALSE.,BAIL)
      END IF
      DO 100 K = 1,4
         ISHIFT(K) = (K-1)*NROWS
  100 CONTINUE
C
      DO 110 I=1,12
         READ(NFT11)
  110 CONTINUE
      READ(NFT11) IARC
      CALL SEQREW(NFT11)
C
      DO 1420 KSTAT = 1,NSTATE
         CALL VCLR(WEIGHT,1,NIRRED)
         NSYMG=0
         IF(NONABEL.AND..NOT.SKIPS) THEN
            CALL VCLR(X(LCCI),1,MAXDET)
            NDET=0
         END IF
         IF (MASWRK) WRITE (IW,9100) KSTAT,EIGVAL(KSTAT)
C
C          THE FOLLOWING CODE SEEMS TO TOUCH 1 EXTRA ELEMENT IN -LEVIR-
C
         IWKS = 0
         LEV = 1
         LEVM = 1
         IR0 = 1
  780    CONTINUE
         IF (LEV .EQ. NLEVS) GO TO 860
         LEVIR(LEV) = IR0
         LEVM = LEV
         LEV = LEVM+1
         LEVIR(LEV) = LEVNR(LEV)+1
  800    CONTINUE
         IR0 = LEVIR(LEV)
         NPTX = LEVPT(LEV)
         NPTM = LEVPT(LEVM)
         IRM = LEVIR(LEVM)
  820    CONTINUE
         IR0 = IR0-1
         IF (IR0 .EQ. 0) GO TO 1000
         NPT = IR0+NPTX
         DO 840 K = 1,4
            IARPT = NPT+ISHIFT(K)
            JARC = IARC(IARPT)
            IF (JARC .EQ. 0) GO TO 840
            JARC = JARC-NPTM
            ICASE(LEVM) = K
            IF (IRM .EQ. JARC) GO TO 780
  840    CONTINUE
         GO TO 820
C
  860    CONTINUE
C
         DO 980 ILEV = 1,NORBS
            ICAS = ICASE(ILEV)
            GO TO (880,900,900,920), ICAS
  880       IOCC = 0
            GO TO 940
C
  900       IOCC = 1
            GO TO 940
C
  920       IOCC = 2
  940       CONTINUE
            IORB = -1
            DO 960 I = 1,NORBMX
               IF (IOUT(I) .EQ. ILEV) IORB = I
  960       CONTINUE
            IF(IORB.GT.0) IECONF(IORB) = IOCC
  980    CONTINUE
C
         IWKS = IWKS+1
         DUM = EIGVEC(IWKS,KSTAT)
         IF((ABS(DUM).GE.PRTTOL) .AND. MASWRK)
     *       WRITE(IW,9120) IWKS,DUM,
     *                  (IECONF(NCORE+IORB),IORB=1,NPART)
C
C        SET SEPARATE TOL FOR SYMMETRY CONSIDERATIONS
C
         IF(RUNTRM.AND.ABELPT()) THEN
           IRREP=IXFTCH(X(LCSFIR),IWKS)
           WEIGHT(IRREP)=WEIGHT(IRREP)+DUM*DUM
         END IF
C
         IF(NONABEL.AND..NOT.SKIPS) THEN
            NONZ=NOIRR.GE.0
            IF(NOIRR.LT.0) THEN
               DO 990 KST=1,NSTATE
C
C                 THIS IS SOMEWHAT QUESTIONABLE PRECISION OF SYMTOL/NT
C                 NOTE THAT EIGVEC GETS MULTIPLIED BY A COEF .LT. ONE
C                 IN CSFTODET TO CONVERT CSF -> DET
C
                  IF(ABS(EIGVEC(IWKS,KST)).GT.SMALLCSF) THEN
                     NONZ=.TRUE.
                     GO TO 995
                  END IF
  990          CONTINUE
  995          CONTINUE
            END IF
            IF(NONZ) CALL CSFTODET(NPART,EIGVEC(IWKS,KSTAT),MXSPIN,NAEL,
     *                             MAXDET,NDET,X(LIGCSF),X(LCCI),ICASE,
     *                             X(LIICAS),X(LIPRIM),X(LCGC),X(LPK),
     *                             X(LNTRAP),X(LBK))
         END IF
 1000    CONTINUE
C
         LEV = LEVM
         LEVM = LEV -1
         IF (LEVM .GT. 0) GO TO 800
C
C           NOW DETERMINE SYMMETRY OF THIS STATE
C
         IF(NONABEL) THEN
            IF(SKIPS) THEN
               IRET=0
            ELSE
               CALL IRRSTATE(IRET,NAO,NCORE,NAEL,MUL,NDET,X(LIGCSF),
     *               X(LCCI),X(LFF),X(LCHARR),X(LIRRID),X(LIRROW),SMALL)
            END IF
            IRREP=IRET
         END IF
C
         IF(RUNTRM) THEN
            IF(ABELPT()) THEN
                IRREP=0
                DO 1050 IG=1,NIRRED
                   IF(WEIGHT(IG).GT.SMALL) THEN
                      NSYMG=NSYMG+1
                      IRREP=IOR(ISHFT(IRREP,4),IG)
                      IF(MASWRK) WRITE(IW,9060) IRPNAM(IG),WEIGHT(IG)
                      IF(MASWRK) CALL FLSHBF(IW)
                   END IF
 1050           CONTINUE
            END IF
            CALL IXSTOR(X(LSTSYM),KSTAT,IRREP)
         END IF
 1420 CONTINUE
C
      IF(RUNTRM) THEN
         DSKWRK=.TRUE.
         CALL RAWRIT(JSODAF,X(JSODA),X(LSTSYM),(NSTATE-1)/NWDVAR+1,
     *              (ICI-1)*NRECJ+6,1)
         DSKWRK=SVDSKW
      END IF
      IF(RUNTRM.OR.NONABEL) CALL RETFM(NEED)
      RETURN
C
 9100 FORMAT(/1X,'STATE #',I5,'  ENERGY =',F18.9//
     *        6X,'CSF',6X,'COEF',4X,'OCCUPANCY (IGNORING CORE)'/
     *        6X,'---',6X,'----',4X,'--------- --------- -----')
 9120 FORMAT(1X,I8,F12.6,2X,50I1,(/23X,50I1))
 9200 FORMAT(/1X,49(1H-)/
     *        1X,'NON-ABELIAN CI WAVEFUNCTION STATE SYMMETRY DRIVER'/
     *        10X,'WRITTEN BY DMITRI FEDOROV.',/,
     *        5X,'MEMORY USED IS',I9,' WORDS'/1X,49(1H-)/)
 9060 FORMAT(/1X,'THE PROJECTION OF THIS CI STATE ONTO SPACE SYMMETRY ',
     *         A4,' WEIGHS',1P,E11.4/)
      END
C*MODULE MTHLIB  *DECK PREV
      SUBROUTINE PREV(V,E,M,N,LDV)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      DIMENSION V(LDV,M),E(M)
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
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
C     ----- PRINT OUT EIGENDATA (VECTORS AND VALUES) -----
C     THE ROWS WILL BE LABELED WITH THE BASIS FUNCTION TAGS.
C     -V- IS N X M, WITH TRUE LEADING DIMENSION -LDV-
C
      IF (MASWRK) THEN
      MAX = 5
      IF (NPRINT .EQ. 6) MAX = 10
      IMAX = 0
C
  100 IMIN = IMAX+1
      IMAX = IMAX+MAX
      IF (IMAX .GT. M) IMAX = M
      WRITE (IW,9008)
      WRITE (IW,9028) (I,I = IMIN,IMAX)
      WRITE (IW,9008)
      WRITE (IW,9068) (E(I),I = IMIN,IMAX)
      WRITE (IW,9008)
      DO 120 J = 1,N
         WRITE (IW,9048) J,BFLAB(J),(V(J,I),I = IMIN,IMAX)
  120 CONTINUE
      IF (IMAX .LT. M) GO TO 100
      END IF
      RETURN
C
 9008 FORMAT(1X)
 9028 FORMAT(15X,10(4X,I4,3X))
 9048 FORMAT(I5,2X,A8,10F11.6)
 9068 FORMAT(15X,10F11.4)
      END
C*MODULE MTHLIB  *DECK PREVNL
      SUBROUTINE PREVNL(V,E,M,N,LDV)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION V(LDV,M),E(M)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- PRINT OUT EIGENDATA (VECTORS AND VALUES) -----
C     THE ROWS ARE NOT LABELED WITH BASIS TAGS.
C     -V- IS N X M, WITH TRUE LEADING DIMENSION -LDV-
C
      IF (MASWRK) THEN
      MAX = 5
      IF (NPRINT .EQ. 6) MAX = 10
      IMAX = 0
  100 IMIN = IMAX+1
      IMAX = IMAX+MAX
      IF (IMAX .GT. M) IMAX = M
      WRITE (IW,9008)
      WRITE (IW,9028) (I,I = IMIN,IMAX)
      WRITE (IW,9008)
      WRITE (IW,9068) (E(I),I = IMIN,IMAX)
      WRITE (IW,9008)
      DO 120 J = 1,N
  120 WRITE (IW,9048) J,(V(J,I),I = IMIN,IMAX)
      IF (IMAX .LT. M) GO TO 100
      END IF
      RETURN
C
 9008 FORMAT(1X)
 9028 FORMAT(15X,10(4X,I4,3X))
 9048 FORMAT(I5,10X,10F11.6)
 9068 FORMAT(15X,10F11.6)
      END
C*MODULE MTHLIB  *DECK PREVS
      SUBROUTINE PREVS(V,E,LABMO,NMO,NAO,LDV,ISTMO)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      DIMENSION V(LDV,NMO),E(NMO),LABMO(NMO)
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
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
C     ----- PRINT OUT EIGENDATA, WITH MO SYMMETRY LABELS -----
C     THE ROWS ARE LABELED WITH THE BASIS FUNCTION NAMES.
C
      IF (MASWRK) THEN
      MAX = 5
      IF (NPRINT .EQ. 6) MAX = 10
      IMAX = ISTMO-1
C
  100 IMIN = IMAX+1
      IMAX = IMAX+MAX
      IF (IMAX .GT. NMO) IMAX = NMO
      WRITE (IW,9008)
      WRITE (IW,9028) (I,       I=IMIN,IMAX)
      WRITE (IW,9068) (E(I),    I=IMIN,IMAX)
      WRITE (IW,9078) (LABMO(I),I=IMIN,IMAX)
      DO 120 J = 1,NAO
         WRITE (IW,9048) J,BFLAB(J),(V(J,I),I = IMIN,IMAX)
  120 CONTINUE
      IF (IMAX .LT. NMO) GO TO 100
      END IF
      RETURN
C
 9008 FORMAT(1X)
 9028 FORMAT(15X,10(4X,I4,3X))
 9048 FORMAT(I5,2X,A8,10F11.6)
 9068 FORMAT(15X,10F11.4)
 9078 FORMAT(16X,10(5X,A4,2X))
      END
C*MODULE MTHLIB  *DECK PRSQ
      SUBROUTINE PRSQ(V,M,N,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION V(NDIM,M)
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
C     ----- PRINT OUT A SQUARE MATRIX -----
C     -V- IS -N- ROWS BY -M- COLUMNS, WITH LEADING DIMENSION -NDIM-
C
      IF (MASWRK) THEN
      MAX = 5
      IF (NPRINT .EQ. 6) MAX = 10
      IMAX = 0
  100 IMIN = IMAX+1
      IMAX = IMAX+MAX
      IF (IMAX .GT. M) IMAX = M
      WRITE (IW,9008)
      WRITE (IW,9028) (I,I = IMIN,IMAX)
      WRITE (IW,9008)
      DO 120 J = 1,N
  120 WRITE (IW,9048) J,(V(J,I),I = IMIN,IMAX)
      IF (IMAX .LT. M) GO TO 100
      END IF
      RETURN
 9008 FORMAT(1X)
 9028 FORMAT(6X,10(4X,I4,4X))
 9048 FORMAT(I5,1X,10F12.7)
      END
C*MODULE MTHLIB  *DECK PRSQL
      SUBROUTINE PRSQL(V,M,N,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      DIMENSION V(NDIM,M)
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
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
C     ----- PRINT OUT A SQUARE MATRIX WITH BASIS SET LABELS -----
C     -V- IS -N- ROWS BY -M- COLUMNS, WITH LEADING DIMENSION -NDIM-
C
      IF (MASWRK) THEN
      MAX = 5
      IF (NPRINT .EQ. 6) MAX = 10
      IMAX = 0
  100 IMIN = IMAX+1
      IMAX = IMAX+MAX
      IF (IMAX .GT. M) IMAX = M
      WRITE (IW,9008)
      WRITE (IW,9028) (I,I = IMIN,IMAX)
      WRITE (IW,9008)
      DO 120 J = 1,N
  120 WRITE (IW,9048) J,BFLAB(J),(V(J,I),I = IMIN,IMAX)
      IF (IMAX .LT. M) GO TO 100
      END IF
      RETURN
 9008 FORMAT(1X)
 9028 FORMAT(15X,10(4X,I4,3X))
 9048 FORMAT(I5,2X,A8,10F11.6)
      END
C*MODULE MTHLIB  *DECK PRTRI
      SUBROUTINE PRTRI(D,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION D(*)
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
C     ----- PRINT SYMMETRIC MATRIX -D- OF DIMENSION -N- -----
C
      IF (MASWRK) THEN
      MAX = 5
      IF (NPRINT .EQ. 6) MAX = 10
      MM1 = MAX-1
      DO 120 I0=1,N,MAX
         IL = MIN(N,I0+MM1)
         WRITE(IW,9008)
         WRITE(IW,9028) (I,I=I0,IL)
         WRITE(IW,9008)
         IL = -1
         DO 100 I=I0,N
            IL=IL+1
            J0=I0+(I*I-I)/2
            JL=J0+MIN(IL,MM1)
            WRITE(IW,9048) I,(D(J),J=J0,JL)
  100    CONTINUE
  120 CONTINUE
      END IF
      RETURN
 9008 FORMAT(1X)
 9028 FORMAT(6X,10(4X,I4,4X))
 9048 FORMAT(I5,1X,10F12.7)
      END
C*MODULE MTHLIB  *DECK PRTRIL
      SUBROUTINE PRTRIL(D,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      DIMENSION D(*)
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
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
C     ----- PRINT SYMMETRIC MATRIX -D- OF DIMENSION -N- -----
C     THE ROWS WILL BE LABELED WITH BASIS FUNCTION TAGS.
C
      IF (MASWRK) THEN
      MAX = 5
      IF (NPRINT .EQ. 6) MAX = 10
      MM1 = MAX - 1
      DO 120 I0=1,N,MAX
         IL = MIN(N,I0+MM1)
         WRITE(IW,9008)
         WRITE(IW,9028) (I,I=I0,IL)
         WRITE(IW,9008)
         IL = -1
         DO 100 I=I0,N
            IL=IL+1
            J0=I0+(I*I-I)/2
            JL=J0+MIN(IL,MM1)
            WRITE(IW,9048) I,BFLAB(I),(D(J),J=J0,JL)
  100    CONTINUE
  120 CONTINUE
      END IF
      RETURN
 9008 FORMAT(1X)
 9028 FORMAT(15X,10(4X,I4,3X))
 9048 FORMAT(I5,2X,A8,10F11.6)
      END
C*MODULE MTHLIB  *DECK PUSQL
      SUBROUTINE PUSQL(V,M,N,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION V(NDIM,M)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- PUNCH A RECTANGULAR MATRIX WITH ORDERING LABELS -----
C     -V- IS -N- ROWS BY -M- COLUMNS, WITH TRUE LEAD DIMENSION -NDIM-
C
      IF (MASWRK) THEN
      DO 120 J = 1,M
      IC = 0
      MAX = 0
  100 MIN = MAX+1
      MAX = MAX+5
      IC = IC+1
      IF (MAX .GT. N) MAX = N
      MODJ=MOD(J,100)
      WRITE (IP,9008) MODJ,IC,(V(I,J),I = MIN,MAX)
      IF (MAX .LT. N) GO TO 100
  120 CONTINUE
      END IF
      RETURN
C
 9008 FORMAT(I2,I3,1P,5E15.8)
      END
C*MODULE MTHLIB  *DECK QMATRX
      SUBROUTINE QMATRX(S,Q,EIG,SCR,IWRK,L0,L1,LDQ,DBG)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION S(*),Q(LDQ,L1),EIG(L1),SCR(L1,8),IWRK(L1)
C
      LOGICAL DBG,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
C
      PARAMETER (ONE=1.0D+00, TOL=1.0D-10,
     *           ZERO=0.0D+00, SMALL=1.0D-11)
C
C     ----- FIND CANONICAL ORTHONORMAL MO-S ------
C     THE "Q-MATRIX" SATISFIES Q-DAGGER * S * Q = I ON EXIT.
C     THE OVERLAP MATRIX -S- IS DESTROYED IN THE DIAGONALIZATION.
C     THE NUMBER OF LINEARLY INDEPENDENT EIGENVECTORS OF -S- IS
C     RETURNED BY -L0-, NOTE THAT THIS ROUTINE WILL KILL THE JOB
C     IN THE CASE THAT -L0- IS NOT EQUAL TO THE INPUT -L1-.
C     THE ARGUMENT -L3- IS UNUSED BY THIS ROUTINE.
C
C        DIAGONALIZE OVERLAP MATRIX
C
      IGERR = 0
      CALL GLDIAG(L1,L1,LDQ,S,SCR,EIG,Q,IGERR,IWRK)
      IF (IGERR .NE. 0) CALL ABRT
C
      IF(DBG) THEN
         WRITE(IW,9000)
         CALL PREVNL(Q,EIG,L1,L1,LDQ)
      END IF
C
C     ----- ELIMINATE EIGENVECTORS WITH SMALL EIGENVALUES -----
C
      DUM = EIG(1)
      J = 0
      K = 0
      DO 180 I = 1,L1
         IF (EIG(I) .LT. TOL) GO TO 160
            J = J+1
            EIG(J) = ONE/SQRT(EIG(I))
            CALL DCOPY(L1,Q(1,I),1,Q(1,J),1)
            GO TO 180
  160    K = K+1
  180 CONTINUE
C
      L0 = L1-K
C
      IF(K.NE.0) THEN
         IF (MASWRK) WRITE(IW,9010) DUM,K,TOL,L0
         CALL ABRT
      END IF
C
C     ----- FORM CANONICAL ORTHONORMAL ORBITALS -----
C
      DO 410 J=1,L0
         EJ = EIG(J)
         DO 400 I=1,L1
            QIJ=Q(I,J)
            QIJ=QIJ*EJ
            IF(ABS(QIJ).LT.SMALL) QIJ=ZERO
            Q(I,J)=QIJ
  400    CONTINUE
  410 CONTINUE
C
      IF(DBG) THEN
         WRITE(IW,9020)
         CALL PREVNL(Q,EIG,L0,L1,LDQ)
      END IF
      RETURN
C
 9000 FORMAT(/5X,17(1H-)/5X,'EIGENVECTORS OF S'/5X,17(1H-))
 9010 FORMAT(20H ..... WARNING ...../
     +     50H THE SMALLEST EIGENVALUE OF THE OVERLAP MATRIX IS ,F15.8/
     +     15H THERE IS(ARE) ,I5,25H EIGENVALUE(S) LESS THAN ,1P,E10.2/
     +     42H THE NUMBER OF CANONICAL ORBITALS KEPT IS ,I5)
 9020 FORMAT(/5X,30(1H-)/5X,'CANONICAL ORTHONORMAL ORBITALS'/5X,30(1H-))
      END
C*MODULE MTHLIB  *DECK ROT
      SUBROUTINE ROT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      COMMON /TRANSF/ XOLD,YOLD,ZOLD,XNEW,YNEW,ZNEW,XP,YP,ZP
      COMMON /FRAME / U1,U2,U3,V1,V2,V3,W1,W2,W3,X0,Y0,Z0
C
C     CALCULATE THE COORDINATES (XP,YP,ZP) OF A POINT IN THE MASTER
C     FRAME GIVEN THE COORDINATES (XNEW,YNEW,ZNEW) IN THE LOCAL FRAME
C
      XP = X0+U1*XNEW+V1*YNEW+W1*ZNEW
      YP = Y0+U2*XNEW+V2*YNEW+W2*ZNEW
      ZP = Z0+U3*XNEW+V3*YNEW+W3*ZNEW
      RETURN
      END
C*MODULE MTHLIB  *DECK RUNITV
      SUBROUTINE RUNITV(NROWS,NCOLS,VECTOR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VECTOR(NROWS,NCOLS)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
C* 23 MAY 1983 - STE
C*
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY - USDOE) MAR 1983
C*
C*    PURPOSE -
C*       GENERATE W.P. REAL UNIT VECTORS
C*                     *    **** *
C*
C*    ON ENTRY -
C*       NROWS  - INTEGER
C*                LENGTH OF VECTOR(S)
C*       NCOLS  - INTEGER
C*                NUMBER OF VECTORS
C*
C*    ON EXIT -
C*       VECTOR - W.P. REAL (NROWS,NCOLS)
C*                UNIT VECTORS ONE TO NCOLS
C*
      DO 120 J=1,NCOLS
         DO 110 I=1,NROWS
            VECTOR(I,J)=ZERO
  110    CONTINUE
         VECTOR(J,J)=ONE
  120 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK SCHMD
      SUBROUTINE SCHMD(V,M,N,LDV)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION V(LDV,N)
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-10)
C
C     ----- ORTHONORMALIZE FIRST -M- MO'S -----
C     AND GENERATE ORTHOGONAL COMPLEMENT SPACE TO FILL
C     OUT A TOTAL OF -N- ORBITALS.
C
      IF (M .EQ. 0) GO TO 180
      DO 160 I = 1,M
      DUMI = ZERO
      DO 100 K = 1,N
  100 DUMI = DUMI+V(K,I)*V(K,I)
      DUMI = ONE/ SQRT(DUMI)
      DO 120 K = 1,N
  120 V(K,I) = V(K,I)*DUMI
C
      IF (I .EQ. M) GO TO 160
      I1 = I+1
      DO 140 J = I1,M
      DUM = -DDOT(N,V(1,J),1,V(1,I),1)
      CALL DAXPY(N,DUM,V(1,I),1,V(1,J),1)
  140 CONTINUE
C
  160 CONTINUE
      IF (M .EQ. N) RETURN
  180 CONTINUE
C
C     ----- GET ORTHOGONAL SPACE -----
C
      I = M
      J = 0
  200 I0 = I
      I = I+1
      IF (I .GT. N) RETURN
  220 J = J+1
      IF (J .GT. N) GO TO 320
      DO 240 K = 1,N
  240 V(K,I) = ZERO
      V(J,I) = ONE
C
      DO 260 II = 1,I0
      DUM = -DDOT(N,V(1,II),1,V(1,I),1)
      CALL DAXPY(N,DUM,V(1,II),1,V(1,I),1)
  260 CONTINUE
      DUMI = ZERO
      DO 280 K = 1,N
  280 DUMI = DUMI+V(K,I)*V(K,I)
      IF ( ABS(DUMI) .LT. TOL) GO TO 220
      DUMI = ONE/ SQRT(DUMI)
      DO 300 K = 1,N
  300 V(K,I) = V(K,I)*DUMI
      GO TO 200
  320 IF (MASWRK) WRITE (IW,9010) I0,N
      CALL ABRT
      STOP
 9010 FORMAT(1X,'ORTHOGONALIZING FOUND A REDUNDANT SET OF VECTORS'/
     *       1X,'ONLY',I5,' INDEPENDENT VECTORS FOUND'/
     *       1X,'BUT THERE SHOULD BE',I5,' TOTAL MOS')
      END
C*MODULE MTHLIB  *DECK SCHMIDT
      SUBROUTINE SCHMIDT(V,S,WRK,N,M,LDV)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION V(LDV,M),S(*),WRK(N)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      PARAMETER (SMALL=1.0D-15)
C
C     GRAM-SCHMIDT ORTHOGONALISATION OF -M- VECTORS -V- OF SIZE -N-
C     ON A METRIC -S- (I.E., OVERLAP INTEGRALS).
C
      NSMALL=0
      DO I=1,M
         CALL MTARBR(S,N,V(1,I),1,WRK,N,1)
         DO J=1,I-1
            CALL DAXPY(N,-DDOT(N,WRK,1,V(1,J),1),V(1,J),1,V(1,I),1)
         ENDDO
         Q=SQRT(DDOT(N,WRK,1,V(1,I),1))
         IF(ABS(Q).LT.SMALL) THEN
            IF(MASWRK) WRITE(IW,9000) I
            NSMALL=NSMALL+1
         ELSE
            CALL DSCAL(N,1.0D+00/Q,V(1,I),1)
         END IF
      ENDDO
      IF(NSMALL.GT.0) CALL ABRT
      RETURN
 9000 FORMAT(1X,'SCHMIDT: VECTOR NORM TOO SMALL FOR ORBITAL',I5)
      END
C*MODULE MTHLIB  *DECK SETHERMITE
      SUBROUTINE SETHERMITE
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /HERMIT/ H11,H21,H22,H31,H32,H33,H41,H42,H43,H44,
     *                H51,H52,H53,H54,H55,H61,H62,H63,H64,H65,H66,
     *                H71,H72,H73,H74,H75,H76,H77,
     *                H81,H82,H83,H84,H85,H86,H87,H88,
     *                H91,H92,H93,H94,H95,H96,H97,H98,H99,
     *                H101,H102,H103,H104,H105,H106,H107,H108,H109,H1010
      COMMON /WERMIT/ W11,W21,W22,W31,W32,W33,W41,W42,W43,W44,
     *                W51,W52,W53,W54,W55,W61,W62,W63,W64,W65,W66,
     *                W71,W72,W73,W74,W75,W76,W77,
     *                W81,W82,W83,W84,W85,W86,W87,W88,
     *                W91,W92,W93,W94,W95,W96,W97,W98,W99,
     *                W101,W102,W103,W104,W105,W106,W107,W108,W109,W1010
C
C     --- INITIALIZE ROOTS AND WEIGHTS FOR GAUSS-HERMITE QUADRATURE ---
C     THE VALUES USED BELOW WERE OBTAINED BY RUNNING THE ROUTINE GIVEN
C     IN THE "NUMERICAL RECIPES" BOOK IN QUADRUPLE PRECISION.  WE WILL
C     LET THE COMPILER TRUNCATE THEM TO THE MACHINE'S DOUBLE PRECISION.
C
      H11=  ZERO
C
      H21= -0.70710678118654752440D+00
      H22=  0.70710678118654752440D+00
C
      H31= -1.22474487139158904910D+00
      H32=  ZERO
      H33=  1.22474487139158904910D+00
C
      H41= -1.65068012388578455588D+00
      H42= -0.52464762327529031788D+00
      H43=  0.52464762327529031788D+00
      H44=  1.65068012388578455588D+00
C
      H51= -2.02018287045608563293D+00
      H52= -0.95857246461381850711D+00
      H53=  ZERO
      H54=  0.95857246461381850711D+00
      H55=  2.02018287045608563293D+00
C
      H61= -2.35060497367449222283D+00
      H62= -1.33584907401369694971D+00
      H63= -0.43607741192761650868D+00
      H64=  0.43607741192761650868D+00
      H65=  1.33584907401369694971D+00
      H66=  2.35060497367449222283D+00
C
      H71= -2.65196135683523349245D+00
      H72= -1.67355162876747144503D+00
      H73= -0.81628788285896466304D+00
      H74= ZERO
      H75=  0.81628788285896466304D+00
      H76=  1.67355162876747144503D+00
      H77=  2.65196135683523349245D+00
C
      H81= -2.93063742025724401922D+00
      H82= -1.98165675669584292585D+00
      H83= -1.15719371244678019472D+00
      H84= -0.38118699020732211685D+00
      H85=  0.38118699020732211685D+00
      H86=  1.15719371244678019472D+00
      H87=  1.98165675669584292585D+00
      H88=  2.93063742025724401922D+00
C
      H91= -3.19099320178152760723D+00
      H91= -2.26658058453184311180D+00
      H91= -1.46855328921666793167D+00
      H91= -0.72355101875283757332D+00
      H91= ZERO
      H91=  0.72355101875283757332D+00
      H91=  1.46855328921666793167D+00
      H91=  2.26658058453184311180D+00
      H91=  3.19099320178152760723D+00
C
      H101=  -3.43615911883773760333D+00
      H102=  -2.53273167423278979641D+00
      H103=  -1.75668364929988177345D+00
      H104=  -1.03661082978951365418D+00
      H105=  -0.34290132722370460879D+00
      H106=   0.34290132722370460879D+00
      H107=   1.03661082978951365418D+00
      H108=   1.75668364929988177345D+00
      H109=   2.53273167423278979641D+00
      H1010=  3.43615911883773760333D+00
C
      W11= 1.77245385090551602730D+00  ! SQRT(PI)
C
      W21= 8.86226925452758013649D-01
      W22= 8.86226925452758013649D-01
C
      W31= 2.95408975150919337883D-01
      W32= 1.18163590060367735153D+00
      W33= 2.95408975150919337883D-01
C
      W41= 8.13128354472451771430D-02
      W42= 8.04914090005512836506D-01
      W43= 8.04914090005512836506D-01
      W44= 8.13128354472451771430D-02
C
      W51= 1.99532420590459132077D-02
      W52= 3.93619323152241159828D-01
      W53= 9.45308720482941881226D-01
      W54= 3.93619323152241159828D-01
      W55= 1.99532420590459132077D-02
C
      W61= 4.53000990550884564086D-03
      W62= 1.57067320322856643916D-01
      W63= 7.24629595224392524092D-01
      W64= 7.24629595224392524092D-01
      W65= 1.57067320322856643916D-01
      W66= 4.53000990550884564086D-03
C
      W71= 9.71781245099519154149D-04
      W72= 5.45155828191270305922D-02
      W73= 4.25607252610127800520D-01
      W74= 8.10264617556807326765D-01
      W75= 4.25607252610127800520D-01
      W76= 5.45155828191270305922D-02
      W77= 9.71781245099519154149D-04
C
      W81= 1.99604072211367619206D-04
      W82= 1.70779830074134754562D-02
      W83= 2.07802325814891879543D-01
      W84= 6.61147012558241291030D-01
      W85= 6.61147012558241291030D-01
      W87= 2.07802325814891879543D-01
      W86= 1.70779830074134754562D-02
      W88= 1.99604072211367619206D-04
C
      W91= 3.96069772632643819046D-05
      W92= 4.94362427553694721722D-03
      W93= 8.84745273943765732880D-02
      W94= 4.32651559002555750200D-01
      W95= 7.20235215606050957124D-01
      W96= 4.32651559002555750200D-01
      W97= 8.84745273943765732880D-02
      W98= 4.94362427553694721722D-03
      W99= 3.96069772632643819046D-05
C
      W101=  7.64043285523262062916D-06
      W102=  1.34364574678123269220D-03
      W103=  3.38743944554810631362D-02
      W104=  2.40138611082314686417D-01
      W105=  6.10862633735325798784D-01
      W106=  6.10862633735325798784D-01
      W107=  2.40138611082314686417D-01
      W108=  3.38743944554810631362D-02
      W109=  1.34364574678123269220D-03
      W1010= 7.64043285523262062916D-06
C
      RETURN
      END
C*MODULE MTHLIB  *DECK SLVLEQ
      SUBROUTINE SLVLEQ(A,X,IPVT,N,LDA,JOB,IERR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,N),X(N),IPVT(N)
C
C     SOLVE LINEAR EQUATIONS.
C     FOR JOB=0, SOLVE AX=B.  FOR JOB.NE.0, SOLVE A-TRANSPOSE*X=B.
C     ON ENTRY, B IS STORED IN X.  ON EXIT, X CONTAINS THE SOLUTION.
C     ON ENTRY, A IS A GENERAL MATRIX.  IT IS DESTROYED ON EXIT.
C     N IS THE SIZE OF THE PROBLEM.  IPVT IS AN INTEGER WORK ARRAY.
C     TWO LINPACK ROUTINES ARE CALLED, TO FACTOR A, THEN SOLVE.
C
      IERR=0
      CALL DGEFA(A,LDA,N,IPVT,IERR)
      IF(IERR.NE.0) RETURN
      CALL DGESL(A,LDA,N,IPVT,X,JOB)
      RETURN
      END
C*MODULE MTHLIB  *DECK SORTIN
      SUBROUTINE SORTIN(IA,LENA)
C
      DIMENSION IA(LENA)
C
C* 28 OCT 1986 - STE * 10 FEB 1984 - STE * 23 MAY 1983 - STE
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY - USDOE)  DATE: 13 MAY 1983
C*
C*    PURPOSE -
C*       SORT AN INTEGER ARRAY  IA  OF LENGTH  LENA  INTO INCREASING
C*       ****    **
C*       ORDER USING A SHELLSORT METHOD.
C*
C*    REFERENCES -
C*       T. BROWN AND G. DIEHR, "SORTING ALGORITHMS FOR MICROCOMPUTERS",
C*       BYTE, P. 482, MAY 1983.
C*
C*    ON ENTRY -
C*       IA     - INTEGER (LENA)
C*                INTEGERS ARRAY TO BE SORTED.
C*       LENA   - INTEGER
C*                NUMBER OF ELEMENTS IN IA TO BE SORTED.
C*
C*    ON EXIT -
C*       IA     - INTEGER (LENA)
C*                ELEMENTS ARE IN INCREASING ORDER.
C*
C*    COMPLEXITY -
C*       LENA**1.26
C
      IF(LENA.LE.1) RETURN
C
C           CALCULATE INITIAL INCREMENT
C
      N=LENA/2
      DO 100 IP2=1,32
         N=N/2
         IF(N.EQ.0) GO TO 110
  100 CONTINUE
  110 CONTINUE
      INC=2**IP2 - 1
C
C              LOOP OVER DECREASING INC
C
      DO 170 N=1,IP2
         IHI=LENA-INC
         DO 160 I=1,IHI
            IF (IA(I) .LE. IA(I+INC)) GO TO 160
C
C                 SWAP ELEMENTS
C
               ITMP=IA(I+INC)
               IA(I+INC) = IA(I)
               IF (I .GT. INC) GO TO 130
                  IA(I) = ITMP
                  GO TO 160
C
  130          CONTINUE
C
C                 CHECK LOWER ELEMENTS
C
               JLO=I-INC
               J = I
               DO 140 JJ=1,JLO,INC
                  J = J - INC
                  IF(ITMP.GE.IA(J)) GO TO 150
                     IA(J+INC) = IA(J)
  140          CONTINUE
               J=J-INC
  150          CONTINUE
               IA(J+INC)=ITMP
C
  160    CONTINUE
         INC=INC/2
  170 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK SYMSQT
      SUBROUTINE SYMSQT(N,A,LDA,T,INCT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LDA,N),T(INCT,N)
      DATA HALF /0.5D+00/
C
C     SYMMETRISE A INTO A TRIANGULAR MATRIX T
C
      IJ=0
      DO 200 I=1,N
         DO 100 J=1,I
            IJ=IJ+1
            T(1,IJ) = (A(I,J)+A(J,I))*HALF
  100    CONTINUE
  200 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK TFDIAG
      SUBROUTINE TFDIAG(H,V,E,WRK,L1,L2,LDV,NMO)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION H(L2),V(LDV,NMO),E(NMO),WRK(L1)
C
      PARAMETER (ZERO=0.0D+00)
C
C     FIND ONLY THE DIAGONAL ELEMENTS -E- OF THE
C     SYMMETRIC MATRIX PRODUCT V-DAGGER * H * V.
C
      DO 400 K=1,NMO
         DO 200 I=1,L1
            DUM = ZERO
            DO 100 J=1,L1
               IJ = (I*I-I)/2 + J
               IF(I.LT.J) IJ = (J*J-J)/2 + I
               DUM = DUM + H(IJ)*V(J,K)
  100       CONTINUE
            WRK(I) = DUM
  200    CONTINUE
         DUM = ZERO
         DO 300 I=1,L1
            DUM = DUM + V(I,K)*WRK(I)
  300    CONTINUE
         E(K) = DUM
  400 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK TFDS
      SUBROUTINE TFDS(DS,NVAR,NCOORD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DS(*)
C
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
C     ----- TRANSFORM THE INTERNAL DISPLACEMENT VECTOR DS TO
C           A CARTESIAN DISPLACEMENT VECTOR  -----
C
      NQ1 = NCOORD
      NQ3 = NQ1*NVAR
C
C     ----- GROW FAST MEMORY -----
C
      LOADFM = 0
      CALL VALFM(LOADFM)
      I10 = 1 + LOADFM
      I20 = I10 + NQ3
      I30 = I20 + NQ1
      LAST = I30
      NEED = LAST - LOADFM
      CALL GETFM(NEED)
C
C     ----- GET THE TRANSFORMATION VECTOR FROM DIRECT ACCESS -----
C
      CALL DAREAD(IDAF,IODA,X(I10),NQ3,50,0)
      DO 120 I = 1,NCOORD
         X(I20-1+I) = DDOT(NVAR,DS,1,X(I10-1+I),NQ1)
  120 CONTINUE
C
C     ----- MOVE THE TRANSFORMED DISPLACEMENT BACK TO DS -----
C
      CALL DCOPY(NCOORD,X(I20),1,DS,1)
C
C     ----- RETURN FAST MEMORY -----
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE MTHLIB  *DECK TFHB
      SUBROUTINE TFHB(HINT,HCART,NVAR,NCOORD,LDI,LDC,NDAF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION HINT(LDI,NVAR),HCART(LDC,NCOORD)
C
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
C     ----- TRANSFORM THE HESSIAN FROM INTERNALS TO CARTESIANS -----
C     NOTE THAT -HINT- AND -HCART- MAY OCCUPY THE SAME STORAGE
C     NVAR   = NUMBER OF INTERNAL COORDS (NVAR.LE.LDI)
C     NCOORD = NUMBER OF CARTESIAN COORDS (NCOORD.LE.LDC)
C     NDAF   = DAF RECORD CONTAINING THE B MATRIX
C
      IF(NVAR.GT.LDI  .OR.  NCOORD.GT.LDC) CALL ABRT
C
      CALL VALFM(LOADFM)
      LB   = 1    + LOADFM
      LHC  = LB   + NVAR*NCOORD
      LHI  = LHC  + (NCOORD*NCOORD+NCOORD)/2
      LWRK = LHI  + (NVAR*NVAR+NVAR)/2
      LAST = LWRK + NVAR
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
C     ----- GET THE B MATRIX TRANSFORMATION -----
C
      CALL DAREAD(IDAF,IODA,X(LB),NVAR*NCOORD,NDAF,0)
C
C     ----- MOVE HESS TO LOWER TRIANGULAR FORM -----
C
      LOC = LHI-1
      DO 110 I = 1,NVAR
         DO 100 J = 1,I
            LOC = LOC+1
            X(LOC) = HINT(I,J)
  100    CONTINUE
  110 CONTINUE
C
C     ----- NOW CALL THE TRANSFORMATION ROUTINE -----
C
      CALL TFTRI(X(LHC),X(LHI),X(LB),X(LWRK),NCOORD,NVAR,NVAR)
C
C     ----- CONVERT INTERNAL HESSIAN TO FULL SQUARE STORAGE -----
C
      LOC = LHC-1
      DO 180 I = 1,NCOORD
         DO 170 J = 1,I
            LOC = LOC + 1
            HCART(I,J) = X(LOC)
            HCART(J,I) = X(LOC)
  170    CONTINUE
  180 CONTINUE
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE MTHLIB  *DECK TFHF
      SUBROUTINE TFHF(HCART,HINT,NVAR,NCOORD,LDI,LDC,NDAF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION HCART(LDC,NCOORD),HINT(LDI,NVAR)
C
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
C     ----- TRANSFORM THE HESSIAN FROM CARTESIANS TO INTERNALS -----
C     NOTE THAT -HINT- AND -HCART- MAY OCCUPY THE SAME STORAGE
C     NVAR   = NUMBER OF INTERNAL COORDS (NVAR.LE.LDI)
C     NCOORD = NUMBER OF CARTESIAN COORDS (NCOORD.LE.LDC)
C
      IF(NVAR.GT.LDI  .OR.  NCOORD.GT.LDC) CALL ABRT
C
      CALL VALFM(LOADFM)
      LBINV = 1     + LOADFM
      LHC   = LBINV + NCOORD*NVAR
      LHI   = LHC   + (NCOORD*NCOORD+NCOORD)/2
      LWRK  = LHI   + (NVAR*NVAR+NVAR)/2
      LAST  = LWRK  + NCOORD
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
C     ----- GET THE B-INVERSE TRANSFORMATION -----
C
      CALL DAREAD(IDAF,IODA,X(LBINV),NCOORD*NVAR,NDAF,0)
C
C     ----- MOVE HESS TO LOWER TRIANGULAR FORM -----
C
      LOC = LHC-1
      DO 110 I = 1,NCOORD
         DO 100 J = 1,I
            LOC = LOC + 1
            X(LOC) = HCART(I,J)
  100    CONTINUE
  110 CONTINUE
C
C     ----- NOW CALL THE TRANSFORMATION ROUTINE -----
C
      CALL TFTRI(X(LHI),X(LHC),X(LBINV),X(LWRK),NVAR,NCOORD,NCOORD)
C
C     ----- CONVERT HESS TO FULL SQUARE STORAGE -----
C
      LOC = LHI-1
      DO 160 I = 1,NVAR
         DO 150 J = 1,I
            LOC = LOC + 1
            HINT(I,J) = X(LOC)
            HINT(J,I) = X(LOC)
  150    CONTINUE
  160 CONTINUE
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE MTHLIB  *DECK TFSQB
      SUBROUTINE TFSQB(V,Q,WRK,M,N,LDQV)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION V(LDQV,M),Q(LDQV,M),WRK(N)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT,PARR
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXSEQ=150, MXCOLS=5)
      PARAMETER (ZERO=0.0D+00, SMALL=1.0D-11)
C
C     ----- BACK TRANSFORM, V = Q*V -----
C     A SET OF ORBITALS IN -Q- IS TRANSFORMED BY ROTATION -V-
C     ON ENTRY, V IS M X M, AND Q IS N X M.  (M.LE.N)
C      ON EXIT, V IS N X M, AND Q IS UNALTERED.
C
C     ----- INITIALIZATION FOR PARALLEL WORK -----
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT  = -1
      L2CNT = -1
      PARR = GOPARR  .AND.  N.GT.MXSEQ
C
      DO 310 J = 1,M,MXCOLS
         JJMAX = MIN(M,J+MXCOLS-1)
C
C     ----- GO PARALLEL! -----
C     TO DECREASE NEXT VALUE OVERHEAD, DO -MXCOLS- COLUMNS AT A TIME
C
         IF(PARR) THEN
            IF (NXT) THEN
               L2CNT = L2CNT + 1
               IF (L2CNT.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
               IF (NEXT.NE.L2CNT) THEN
                  DO 010 JJ=J,JJMAX
                     CALL VCLR(V(1,JJ),1,N)
  010             CONTINUE
                  GO TO 310
               END IF
            ELSE
               IPCOUNT = IPCOUNT + 1
               IF (MOD(IPCOUNT,NPROC).NE.0) THEN
                  DO 020 JJ=J,JJMAX
                     CALL VCLR(V(1,JJ),1,N)
  020             CONTINUE
                  GO TO 310
               END IF
            END IF
         END IF
C
         DO 300 JJ=J,JJMAX
            DO 100 I = 1,N
               W = DDOT(M,Q(I,1),LDQV,V(1,JJ),1)
               IF(ABS(W).LT.SMALL) W=ZERO
               WRK(I)=W
  100       CONTINUE
            DO 200 I = 1,N
               V(I,JJ) = WRK(I)
  200       CONTINUE
  300    CONTINUE
  310 CONTINUE
C
      IF(PARR) THEN
         IF(LDQV.GT.N) THEN
            NP1 = N + 1
            DO 420 I=NP1,LDQV
               DO 410 J=1,M
                  V(I,J) = ZERO
  410          CONTINUE
  420       CONTINUE
         END IF
         CALL DDI_GSUMF(510,V,LDQV*M)
         IF(NXT) CALL DDI_DLBRESET
      END IF
      RETURN
      END
C*MODULE MTHLIB  *DECK TFSQU
      SUBROUTINE TFSQU(H,F,T,WRK,N,M)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION H(M,M),F(N,N),T(N,M),WRK(N)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT,PARR
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXSEQ=150, MXROWS=5)
C
C     ----- TRANSFORM THE SQUARE MATRIX F USING VECTORS T -----
C                      H = T-DAGGER * F * T
C     THE ORDER OF THE SQUARE MATRICES H, F, AND T IS N.
C
C     ----- INITIALIZATION FOR PARALLEL WORK -----
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT  = -1
      L2CNT = -1
      PARR = GOPARR  .AND.  N.GT.MXSEQ
C
      DO 310 I = 1,M,MXROWS
         IIMAX = MIN(M,I+MXROWS-1)
C
C     ----- GO PARALLEL! -----
C     TO DECREASE NEXT VALUE OVERHEAD, WE DO -MXROWS- AT A TIME
C
         IF(PARR) THEN
            IF (NXT) THEN
               L2CNT = L2CNT + 1
               IF (L2CNT.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
               IF (NEXT.NE.L2CNT) THEN
                  DO 010 II=I,IIMAX
                     CALL VCLR(H(II,1),M,M)
  010             CONTINUE
                  GO TO 310
               END IF
            ELSE
               IPCOUNT = IPCOUNT + 1
               IF (MOD(IPCOUNT,NPROC).NE.0) THEN
                  DO 020 II=I,IIMAX
                     CALL VCLR(H(II,1),M,M)
  020             CONTINUE
                  GO TO 310
               END IF
            END IF
         END IF
C
         DO 300 II=I,IIMAX
            DO 100 L=1,N
               WRK(L) = DDOT(N,T(1,II),1,F(1,L),1)
  100       CONTINUE
            DO 200 J=1,M
               H(II,J) = DDOT(N,WRK,1,T(1,J),1)
  200       CONTINUE
  300    CONTINUE
  310 CONTINUE
C
      IF(PARR) THEN
         CALL DDI_GSUMF(515,H,M*M)
         IF(NXT) CALL DDI_DLBRESET
      END IF
C
      RETURN
      END
C*MODULE MTHLIB  *DECK TFTRI
      SUBROUTINE TFTRI(H,F,T,WRK,M,N,LDT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION H(*),F(*),T(LDT,M),WRK(N)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT,PARR
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXSEQ=150, MXROWS=5)
      PARAMETER (ZERO=0.0D+00, SMALL=1.0D-11)
C
C     ----- TRANSFORM THE TRIANGULAR MATRIX F USING VECTORS T -----
C                      H = T-DAGGER * F * T
C     THE ORDER OF THE TRIANGULAR MATRICES H AND F ARE M AND N.
C
C     ----- INITIALIZATION FOR PARALLEL WORK -----
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT  = -1
      L2CNT = -1
      PARR = GOPARR  .AND.  N.GT.MXSEQ
C
      IF(PARR) THEN
         M2 = (M*M+M)/2
         CALL VCLR(H,1,M2)
      END IF
C
      IJ = 0
      DO 310 J = 1,M,MXROWS
         JJMAX = MIN(M,J+MXROWS-1)
C
C     ----- GO PARALLEL! -----
C
         IF(PARR) THEN
            IF (NXT) THEN
               L2CNT = L2CNT + 1
               IF (L2CNT.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
               IF (NEXT.NE.L2CNT) THEN
                  DO 010 JJ=J,JJMAX
                     IJ = IJ + JJ
  010             CONTINUE
                  GO TO 310
               END IF
            ELSE
               IPCOUNT = IPCOUNT + 1
               IF (MOD(IPCOUNT,NPROC).NE.0) THEN
                  DO 020 JJ=J,JJMAX
                     IJ = IJ + JJ
  020             CONTINUE
                  GO TO 310
               END IF
            END IF
         END IF
C
C             FIRST CALCULATE T-DAGGER TIMES -F-, A ROW AT A TIME
C
         DO 300 JJ=J,JJMAX
            IK = 0
            DO 140 I = 1,N
               IM1 = I-1
               DUM = ZERO
               TDUM = T(I,JJ)
               IF (IM1.GT.0) THEN
                  DO 100 K = 1,IM1
                     IK = IK+1
                     WRK(K) = WRK(K)+F(IK)*TDUM
                     DUM = DUM+F(IK)*T(K,JJ)
  100             CONTINUE
               END IF
               IK = IK+1
               WRK(I) = DUM+F(IK)*TDUM
  140       CONTINUE
C
C             THEN TAKE THAT ROW TIMES EVERY COLUMN IN -T-
C
            DO 200 I = 1,JJ
               IJ = IJ+1
               HIJ = DDOT(N,T(1,I),1,WRK,1)
               IF(ABS(HIJ).LT.SMALL) HIJ=ZERO
               H(IJ)=HIJ
  200       CONTINUE
  300    CONTINUE
  310 CONTINUE
C
      IF(PARR) THEN
         CALL DDI_GSUMF(520,H,M2)
         IF(NXT) CALL DDI_DLBRESET
      END IF
C
      RETURN
      END
C*MODULE MTHLIB  *DECK TFTRIB
      SUBROUTINE TFTRIB(FAO,FMO,S,V,SV,WRK,L0,L1,L2,L3)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION FAO(L2),FMO(L2),S(L2),V(L3),SV(L3),WRK(L1)
C
C        BACK-TRANSFORM A SYMMETRIC OPERATOR EXPRESSED IN
C        THE MO BASIS -V- TO THE AO BASIS, ACCORDING TO
C           FAO = SV * FMO * (SV)-DAGGER
C        THIS IS THE INVERSE TO THE -TFTRI- ROUTINE.
C
C        -S-, -V-, AND -FMO- MUST BE PROVIDED ON ENTRY, AND ARE
C        UNCHANGED ON EXIT.  -SV- AND -WRK- ARE WORKING STORAGE.
C        FOR THE PRESENT, ALL SIZES ARE THE SAME, -V- IS SQUARE.
C
      CALL MTARBR(S,L1,V,L0,SV,L1,1)
      CALL TRPOSQ(SV,L1)
      CALL TFTRI(FAO,FMO,SV,WRK,L1,L0,L1)
      RETURN
      END
C*MODULE MTHLIB  *DECK TRACEP
      DOUBLE PRECISION FUNCTION TRACEP(A,B,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(*),B(*)
C
C     ----- COMPUTE 1 ELECTRON EXPECTATION VALUE ----
C     THE DENSITY MATRIX IS STORED IN A IN TRIANGULAR ORDER.
C     NOTE THAT THE OFF DIAGONAL ELEMENTS OF THE DENSITY ARE
C     TO BE DOUBLED.
C     THE PROPERTY INTEGRALS ARE IN B IN TRIANGULAR ORDER.
C     THE EXPECTATION VALUE IS THE DOT PRODUCT OF THESE MATRICES.
C
      N2 = (N*N+N)/2
      TRACE = DDOT(N2,A,1,B,1)
C                 DOUBLE EVERYTHING, THEN SUBTRACT THE DIAGONAL
      TRACE = TRACE+TRACE
      K = 0
      DO 100 I = 1,N
         K = K+I
         TRACE = TRACE - A(K)*B(K)
  100 CONTINUE
      TRACEP = TRACE
      RETURN
      END
C*MODULE MTHLIB  *DECK TRANG
      SUBROUTINE TRANG(G,NVAR,NCOORD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION G(*)
C
      LOGICAL OUT,GOPARR,DSKWRK,MASWRK
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
      COMMON /FMCOM / X(1)
C
C     ----- TRANSFORM THE GRADIENT FROM CARTESIANS TO INTERNALS -----
C     G(INT) = BINV-DAGGER * G(CART)
C
      OUT = NPRINT.EQ.1 .AND. MASWRK
      NQ1 = NCOORD
      NQ3 = NQ1*NVAR
C
C     ----- GROW FAST MEMORY -----
C
      CALL VALFM(LOADFM)
      I10  = 1 + LOADFM
      I20  = I10 + NQ3
      LAST = I20 + NVAR
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
C     ----- TRANSFORM TO INTERNALS USING BINV MATRIX -----
C
      CALL DAREAD(IDAF,IODA,X(I10),NQ3,50,0)
      IADD = I10 - NQ1
      DO 100 I = 1,NVAR
         IADD = IADD + NQ1
         X(I20-1+I) = DDOT(NCOORD,G,1,X(IADD),1)
  100 CONTINUE
C
C     ----- MOVE THE TRANSFORMED GRADIENT BACK INTO G -----
C
      CALL DCOPY(NVAR,X(I20),1,G,1)
      IF(OUT) WRITE(IW,9008) (G(I),I=1,NVAR)
      CALL RETFM(NEED)
      RETURN
C
 9008 FORMAT(1X,'TRANG: TRANSFORMED GRADIENT'/(6E14.7/))
      END
C*MODULE MTHLIB  *DECK TRANGB
      SUBROUTINE TRANGB(G,NVAR,NCOORD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION G(*)
C
      LOGICAL OUT,GOPARR,DSKWRK,MASWRK
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
      COMMON /FMCOM / X(1)
C
C     ----- TRANSFORM THE GRADIENT FROM INTERNALS TO CARTESIANS -----
C
      OUT = NPRINT .EQ. 1 .AND. MASWRK
      NQ1 = NVAR
      NQ3 = NQ1*NCOORD
C
C     ----- GROW FAST MEMORY -----
C
      LOADFM = 0
      CALL VALFM(LOADFM)
      I10 = 1 + LOADFM
      I20 = I10 + NQ3
      I30 = I20 + NCOORD
      LAST = I30
      NEED = LAST - LOADFM
      CALL GETFM(NEED)
C
C     ----- GET THE TRANSFORMATION VECTORS -----
C
      CALL DAREAD(IDAF,IODA,X(I10),NQ3,49,0)
      IADD = I10 - NQ1
      DO 100 I = 1,NCOORD
         IADD = IADD + NQ1
         X(I20-1+I) = DDOT(NVAR,G,1,X(IADD),1)
  100 CONTINUE
C
C     ----- MOVE THE TRANSFORMED GRADIENT BACK INTO G -----
C
      CALL DCOPY(NCOORD,X(I20),1,G,1)
      CALL DAWRIT(IDAF,IODA,X(I20),NCOORD,3,0)
      IF (OUT) WRITE (IW,9008) (X(I20-1+I),I = 1,NCOORD)
C
C     ----- RETURN FAST MEMORY -----
C
      CALL RETFM(NEED)
      RETURN
 9008 FORMAT(19H PURIFIED GRADIENT ,/,10(6E14.7,/))
      END
C*MODULE MTHLIB  *DECK TRANS
      SUBROUTINE TRANS(NN)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000,MXSH=5000)
C
      COMMON /TRANSF/ XOLD,YOLD,ZOLD,XNEW,YNEW,ZNEW,XP,YP,ZP
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
C      CALCULATE THE COORDINATES (XNEW,YNEW,ZNEW) OF THE TRANSFORM
C      OF THE POINT (XOLD,YOLD,ZOLD) UNDER THE TRANSFORMATION T
C
      XNEW = XOLD*T(NN+1)+YOLD*T(NN+2)+ZOLD*T(NN+3)
      YNEW = XOLD*T(NN+4)+YOLD*T(NN+5)+ZOLD*T(NN+6)
      ZNEW = XOLD*T(NN+7)+YOLD*T(NN+8)+ZOLD*T(NN+9)
      RETURN
      END
C*MODULE MTHLIB  *DECK TRPOSE
      SUBROUTINE TRPOSE(A,B,N,M,KIND)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(N,M), B(M,N)
C
C* 14 JAN 1983 - STE * 8 MAR 1980
C*
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY-USDOE)
C*
C*    PURPOSE -
C*       STORE TRANSPOSE OF N BY M MATRIX A IN MATRIX B OR A
C*             **   ****
C*
C*    ON ENTRY -
C*       A     - W.P. REAL (N,M)
C*               MATRIX TO BE TRANSPOSED
C*       N      - INTEGER
C*                ROWS OF INPUT MATRIX, COLUMNS OF OUTPUT MATRIX
C*       M      - INTEGER
C*                COLUMNS OF INPUT MATRIX, ROWS OF OUTPUT MATRIX
C*       KIND   - INTEGER
C*                IF NOT ZERO, TRANSPOSED MATRIX IS COPYIED BACK INTO A
C*
C*    ON EXIT -
C*       B      - W.P. REAL (M,N)
C*                TRANSPOSED COPY OF INPUT MATRIX
C*       A (OPTIONAL) - W.P. REAL (M,N)
C*                TRANSPOSED COPY OF INPUT MATRIX
C
      IF(N.LE.0 .OR. M.LE.0) RETURN
      DO 120 J=1,M
         DO 110 I=1,N
            B(J,I) = A(I,J)
  110    CONTINUE
  120 CONTINUE
      IF(KIND.NE.0) CALL DCOPY(M*N,B,1,A,1)
      RETURN
      END
C*MODULE MTHLIB  *DECK TRPOSQ
      SUBROUTINE TRPOSQ(A,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(N,N)
C
C     TRANSPOSE SQUARE MATRIX IN PLACE
C
      DO 100 J = 2,N
         JMO = J - 1
         DO 100 I = 1,JMO
            TMP = A(I,J)
            A(I,J) = A(J,I)
            A(J,I) = TMP
  100 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK ZEROD
      SUBROUTINE ZEROD(V,D,E,L1,L2,L3)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION E(L1),D(L2),V(L3)
C
C     ----- ZERO THE WAVEFUNCTION AND ITS DENSITY -----
C
      CALL VCLR(E,1,L1)
      CALL VCLR(D,1,L2)
      CALL VCLR(V,1,L3)
      RETURN
      END
C*MODULE MTHLIB  *DECK VADD
      SUBROUTINE VADD(A,I,B,J,C,K,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(*),B(*),C(*)
C
C     ----- COMPUTE VECTOR C = A+B -----
C     I,J,K ARE INCREMENTS THROUGH A,B,C RESPECTIVELY
C
      IF (I .NE. 1) GO TO 200
      IF (J .NE. 1) GO TO 200
      IF (K .NE. 1) GO TO 200
C
C        I = J = K = 1
C
      DO 110 L=1,N
         C(L) = A(L) + B(L)
  110 CONTINUE
      RETURN
C
  200 CONTINUE
      LA=1-I
      LB=1-J
      LC=1-K
      DO 210 L=1,N
         LA=LA+I
         LB=LB+J
         LC=LC+K
         C(LC) = A(LA) + B(LB)
  210 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK VCLR
      SUBROUTINE VCLR(A,INCA,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(*)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- ZERO OUT VECTOR -A-, USING INCREMENT -INCA- -----
C
      IF (INCA .NE. 1) GO TO 200
      DO 110 L=1,N
         A(L) = ZERO
  110 CONTINUE
      RETURN
C
  200 CONTINUE
      LA=1-INCA
      DO 210 L=1,N
         LA=LA+INCA
         A(LA) = ZERO
  210 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK VSUB
      SUBROUTINE VSUB(A,I,B,J,C,K,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(*),B(*),C(*)
C
C     ----- COMPUTE VECTOR C = B-A -----
C     I,J,K ARE INCREMENTS THROUGH A,B,C RESPECTIVELY
C
      IF (I .NE. 1) GO TO 200
      IF (J .NE. 1) GO TO 200
      IF (K .NE. 1) GO TO 200
C
C        I = J = K = 1
C
      DO 110 L=1,N
         C(L) = B(L) - A(L)
  110 CONTINUE
      RETURN
C
  200 CONTINUE
      LA=1-I
      LB=1-J
      LC=1-K
      DO 210 L=1,N
         LA=LA+I
         LB=LB+J
         LC=LC+K
         C(LC) = B(LA) - A(LB)
  210 CONTINUE
      RETURN
      END
C
C*MODULE MTHLIB  *DECK VICLR
      SUBROUTINE VICLR(IA,INCA,N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION IA(*)
C
C     ----- ZERO OUT AN INTEGER VECTOR -A-, USING INCREMENT -INCA- -----
C
      IF (INCA .NE. 1) GO TO 200
      DO 110 L=1,N
         IA(L) = 0
  110 CONTINUE
      RETURN
C
  200 CONTINUE
      LA=1-INCA
      DO 210 L=1,N
         LA=LA+INCA
         IA(LA) = 0
  210 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK CLEBSCH
      FUNCTION CLEBSCH(J1,J2,M1,M2,J,M)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      PARAMETER(ZERO=0.0D+00,ONE=1.0D+00)
      COMMON /FACTS/FACT(0:100),NFACT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     CALCULATE THE CLEBSCH-GORDAN COEFFICIENTS FOR COUPLING OF
C     TWO ARBITRARY ANGULAR MOMENTA, USING RACAH FORMULA FROM:
C     "ANGULAR MOMENTUM", D.M.BRINK,G.R.SATCHLER,OXFORD,1968
C     <J1 J2 M1 M2 | J M>
C     INPUT: -DOUBLED- VALUES OF ALL ARGUMENTS AS INTEGERS:
C     I.E.: 2*J1, 2*J2, 2*M1, 2*M2, 2*J, 2*M
C     WRITTEN BY DMITRI FEDOROV
C
      CLEBSCH=ZERO
      IF(ABS(M1).GT.J1.OR.ABS(M2).GT.J2.OR.ABS(M).GT.J.OR.J1.LT.0.OR.
     *   J2.LT.0.OR.J.LT.0.OR.ABS(J1-J2).GT.J.OR.J.GT.J1+J2.OR.M1+M2.
     *   NE.M) RETURN
      IF(MAX(J1,J2)+J.GT.NFACT) THEN
         IF(MASWRK) THEN
         WRITE(IW,*) 'INCREASE THE NUMBER OF FACTORIALS PRECALCULATED'
         WRITE(IW,*) 'HAVE ',NFACT,' NEED ',MAX(J1,J2)+J
         END IF
         CALL ABRT
      END IF
      J1NM1=(J1-M1)/2
      JNJ2PM1=(J-J2+M1)/2
      J2PM2=(J2+M2)/2
      JNJ1NM2=(J-J1-M2)/2
      J1PJ2NJ=(J1+J2-J)/2
C     CHECK IF THE EVENNESS IS VALID (EG J1 AND M1 BE BOTH EVEN OR ODD)
      IF(J1NM1*2.NE.J1-M1.OR.J2PM2*2.NE.J2+M2.OR.J1PJ2NJ*2.NE.J1+J2-J)
     * RETURN
      MINK=MAX(MAX(-JNJ2PM1,-JNJ1NM2),0)
      MAXK=MIN(MIN(J1NM1,J2PM2),J1PJ2NJ)
      IF((MINK/2)*2.NE.MINK) THEN
         IPHASE=-1
      ELSE
         IPHASE=1
      END IF
      DO 100 K=MINK,MAXK
         CLEBSCH=CLEBSCH+IPHASE/(FACT(J1NM1-K)*FACT(JNJ2PM1+K)*
     *           FACT(J2PM2-K)*FACT(JNJ1NM2+K)*FACT(K)*FACT(J1PJ2NJ-K))
         IPHASE=-IPHASE
  100 CONTINUE
      IF(MINK.GT.MAXK) CLEBSCH=ONE
      CLEBSCH=CLEBSCH*SQRT(FACT(J1PJ2NJ)*FACT((J1+J-J2)/2)*
     *        FACT((J2+J-J1)/2)/FACT((J1+J2+J)/2+1)*(J+1)*
     *        FACT((J1+M1)/2)*FACT(J1NM1)*FACT(J2PM2)*FACT((J2-M2)/2)*
     *        FACT((J+M)/2)*FACT((J-M)/2))
      RETURN
      END
C*MODULE MTHLIB  *DECK IEULER
      SUBROUTINE IEULER(IEULERR,A,ALPHA,BETA,GAMMA,COSB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      DIMENSION A(3,3)
C     EQUIVALENCE (ANGLES(1),ALPHA),(ANGLES(2),BETA),(ANGLES(3),GAMMA),
C    *            (ANGLES(4),COSB)
      PARAMETER(ZERO=0.0D+00,ONE=1.0D+00,PI=3.14159265358979312D+00,
     *          TWO=2.0D+00)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     INPUT: ANY 3X3 ORTHOGONAL MATRIX
C     OUTPUT: THREE EULER ANGLES AND PARITY 1 OR -1
C     FOR THE DOUBLE GROUPS THIS ROUTINE DOES NOT GENERATE THE STANDARD
C     ANGLES (IE. C4Z**-1 SHOULD BE 4PI-2PI/4, NOT 2PI-2PI/4).
C     THE ANGLES NEED NOT BUT WILL BE FIXED IN ROUTINE KINE.
C
      TOL=1.0D-12
C     SEE TOL IN ROUTINE SAMBO
      DETA=DET3(A)
      IF(ABS(ABS(DETA)-ONE).GT.TOL) THEN
         IF(MASWRK) THEN
         WRITE(IW,*) 'NON-ORTHOGONAL MATRIX IN EULER.'
         CALL PRSQ(A,3,3,3)
         END IF
         CALL ABRT
      END IF
      IEULERR=1
      IF(ABS(DETA-ONE).GT.TOL) IEULERR=-1
      PARITY=IEULERR
      COSB=A(3,3)*PARITY
      BETA=ACOS1(COSB)
      SINB=SIN(BETA)
      IF(ABS(SINB).LE.TOL) THEN
C        AMBIGUITY IN SPLITTING ALPHA+GAMMA
C        VOLUNTARILY SET GAMMA TO 0 AND ASSIGN ALPHA+GAMMA TO ALPHA
         GAMMA=ZERO
         ALPHA=ACOS1(COSB*A(1,1)*PARITY)
C        IF(A(2,1)*COSB*PARITY.LT.ZERO) ALPHA=TWO*PI-ALPHA
         IF(A(2,1)*COSB*PARITY.LT.-TOL) ALPHA=TWO*PI-ALPHA
      ELSE
         GAMMA=ACOS1(-A(3,1)*PARITY/SINB)
C        IF(A(3,2)*PARITY.LT.ZERO) GAMMA=TWO*PI-GAMMA
         IF(A(3,2)*PARITY.LT.-TOL) GAMMA=TWO*PI-GAMMA
         ALPHA=ACOS1(A(1,3)*PARITY/SINB)
C        IF(A(2,3)*PARITY.LT.ZERO) ALPHA=TWO*PI-ALPHA
         IF(A(2,3)*PARITY.LT.-TOL) ALPHA=TWO*PI-ALPHA
      END IF
      RETURN
      END
C*MODULE MTHLIB  *DECK INIFAC
      SUBROUTINE INIFAC(N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      PARAMETER(ONE=1.0D+00)
      COMMON /FACTS/FACT(0:100),NFACT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     FILL THE ARRAY WITH FACTORIALS
      IF(N.GT.100) THEN
         IF(MASWRK) THEN
         WRITE(IW,*) 'CANNOT CALCULATE MORE THAN 100 FACTORIALS.'
         WRITE(IW,*) 'PLEASE INCREASE THE SIZE OF THE INTERNAL ARRAY.'
         END IF
         CALL ABRT
      END IF
      FACT(0)=ONE
      DO 100 I=1,N
         FACT(I)=FACT(I-1)*I
  100 CONTINUE
      NFACT=N
      RETURN
      END
C*MODULE MTHLIB  *DECK ACOS1
      FUNCTION ACOS1(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(ZERO=0.0D+00,ONE=1.0D+00,PI=3.14159265358979312D+00)
C
C     SAME AS ACOS(X) BUT ALLOWS FOR ARGUMENTS OF 1+TOL AND -1-TOL
      TOL=1.0D-12
      IF(ABS(X-ONE).LT.TOL) THEN
         ACOS1=ZERO
      ELSE IF(ABS(X+ONE).LT.TOL) THEN
         ACOS1=PI
           ELSE
              ACOS1=ACOS(X)
           END IF
      RETURN
      END
C*MODULE MTHLIB  *DECK DET3
      FUNCTION DET3(A)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(3,3)
C
C     DETERMINANT OF A 3X3 MATRIX (WHAT ELSE?)
C
      DET3 = A(1,1)*A(2,2)*A(3,3) + A(1,2)*A(2,3)*A(3,1)
     *     + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)
     *     - A(1,1)*A(2,3)*A(3,2) - A(1,2)*A(2,1)*A(3,3)
      RETURN
      END
C*MODULE MTHLIB  *DECK DACOPY
      SUBROUTINE DACOPY(N,DA,DX,INCX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION DX(*)
C
C     INITIALISES A VECTOR WITH A CONSTANT.
C           DX(I) <== DA
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C     CLONED FROM DCOPY BY D. FEDOROV
C
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      DO 10 I = 1,N
        DX(IX) = DA
        IX = IX + INCX
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,7)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        DX(I)=DA
   30 CONTINUE
      IF( N .LT. 7 ) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
        DX(I)=DA
        DX(I + 1)=DA
        DX(I + 2)=DA
        DX(I + 3)=DA
        DX(I + 4)=DA
        DX(I + 5)=DA
        DX(I + 6)=DA
   50 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK SYMTRZE
      SUBROUTINE SYMTRZE(A,LDA,N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,N)
      PARAMETER (HALF=0.5D+00)
C
C     ---- SYMMETRIZE A SQUARE MATRIX A (N X N) ----
C
      DO 200 I=1,N
         DO 100 J=1,I
            VAL = (A(I,J)+A(J,I))*HALF
            A(I,J)=VAL
            A(J,I)=VAL
  100    CONTINUE
  200 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK ASYMTRZE
      SUBROUTINE ASYMTRZE(A,LDA,N,ASYM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,N)
      PARAMETER (HALF=0.5D+00,ZERO=0.0D+00)
C
C     ---- ANTISYMMETRIZE A SQUARE MATRIX A (N X N) ----
C     AND CALCULATE THE AVERAGE A-ANTISYMMMETRY
C     (HOW MUCH NON-ANTISYMMETRIC)
C
      ASYM=ZERO
      DO 200 I=1,N
         DO 100 J=1,I
            VAL = (A(I,J)-A(J,I))*HALF
            ASYM= ASYM + ABS(A(I,J)+A(J,I))
            A(I,J)=VAL
            A(J,I)=-VAL
  100    CONTINUE
  200 CONTINUE
      N2 = (N*N+N)/2
      ASYM = ASYM/N2
C
C     INCLUDING THE DIAGONAL (THAT SHOULD BE ZERO)
C
      RETURN
      END
C*MODULE MTHLIB  *DECK IXFTCH
      INTEGER FUNCTION IXFTCH(IA,I)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION IA(*)
C
C         FETCH AN INTEGER VALUE LURKING IN FLOATING POINT STORAGE,
C         FROM THE POINT OF VIEW OF THE CALLING PROGRAM.
C
      IXFTCH=IA(I)
      RETURN
      END
C*MODULE MTHLIB  *DECK IXSTOR
      SUBROUTINE IXSTOR(IA,I,IELEM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION IA(*)
C
C         STORE AN INTEGER VALUE INTO AN ARRAY.  THIS CALL IS USEFUL
C         IF THE CALLER'S -IA- ARRAY IS DYNAMIC MEMORY, WHICH APPEARS
C         TO BE OF FLOATING POINT TYPE, PREVENTING AN ASSIGNMENT
C         STATEMENT TO BE USED.  SEE -IXFTCH- FOR THE INVERSE OPERATION.
C
      IA(I)=IELEM
      RETURN
      END
C*MODULE MTHLIB  *DECK IRRSTATE
      SUBROUTINE IRRSTATE(IRRSTAT,NAO,NCORE,NAEL,MUL,MAXDET,IDTCSF,
     *                    CDET,FF,CHARR,IRRID,IRROW,SMALL)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,DEBUG
      PARAMETER (MXATM=2000,MXSH=5000,MAXL=3,ZERO=0.0D+00,ONE=1.0D+00)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,NSAFMO
      COMMON /SYMREP/ IRPNAM(14),IPA(14),LAMBDA(14),LAMBD0(14),
     *                IADDR1(14),IADDR2(14),IADDR3(14)
      COMMON /GUGWFN/ NFZC,NMCC,NDOC,NAOS,NBOS,NALP,NVAL,NEXT,NFZV,
     *                IFORS,IEXCIT,ICICI,NOIRR
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      DIMENSION IDTCSF(NAEL,MAXDET),CDET(MAXDET),
     *          FF(MAXL,NT,NAO),CHARR(NT,NIRRED),IRRID(NAO),IRROW(NAO)
      DATA DBGME/8HIRRSTATE/,DEBUGG/8HDEBUG   /
C
C     THIS ROUTINE DETERMINES THE STATE SYMMETRY OF CI WAVEFUNCTIONS FOR
C     NON-ABELIAN GROUPS
C
      DEBUG=(EXETYP.EQ.DBGME.OR.EXETYP.EQ.DEBUGG).AND.MASWRK
      IRRSTAT=0
      NSYMG=0
      TINY=SMALL*SMALL
C     TOL=1.0D-04
C     TOLPRT=5.0D-02
      KCRIT=(NAEL+MUL-1)/2
C
C     THIS IS THE LAST NUMBER OF AN ALPHA ORBITAL IN A DETERMINANT,
C     ASSUMING MS=S!!
C
      IF(DEBUG) THEN
         WRITE(IW,*) 'DECOMPOSED INTO DETERMINANTS:'
         DO 400 I=1,MAXDET
            WRITE(IW,9000) I,CDET(I),(IDTCSF(J,I),J=1,NAEL)
  400    CONTINUE
      END IF
      TOTWEIGHT=ZERO
      DO 900 IG=1,NIRRED
        WEIGHT=ZERO
        DO 800 I=1,MAXDET
          IF(MOD(I,NPROC).NE.ME) GO TO 800
C
C         DETERMINANTS IN LOOP OVER I ARE BUILT UPON SMO(S) (ROW
C         SYMMETRY-ADAPTED MOS), THESE ARE OBTAINED UPON APPLICATION
C         OF PROJECTION OPERATOR ON DETERMINANTS IN THE LOOP OVER J
C         BUILT UPON MOS.
C         IT IS ASSUMED THAT SMO-DETERMINANTS ARE EXACTLY THE SAME IN
C         TERMS OF ORBITAL OCCUPATIONS AS MO-DETERMINANTS.
C         IT MUST BE SO FOR SANE ACTIVE SPACES.
C
          SUM=ZERO
          DO 700 J=1,MAXDET
            IF(IDETDISCO(NCORE,NAEL,IDTCSF(1,I),IDTCSF(1,J),IRRID)
     *            .EQ.0) GO TO 700
            IF(DEBUG) WRITE(IW,*) 'COMPARING DETS',I,J
            SUM1=ZERO
            DO 600 IT=1,NT
              FACT=CHARR(IT,IG)
C             CONJUGATED
              K=1
  500         CONTINUE
C             DO 500 K=1,NAEL
                IMOR=IRROW(ABS(IDTCSF(K,I))-NCORE)
                JMO=ABS(IDTCSF(K,J))-NCORE
                IF(K.LT.NAEL) THEN
C                 IMORNEXT=IRROW(ABS(IDTCSF(K+1,I))-NCORE)
                  JMONEXT=ABS(IDTCSF(K+1,J))-NCORE
C                 THIS TAKES CARE OF THE ONLY CASE OF PHASE APPEARANCE
C                 FOR DOUBLY-DEGENERATE ORBITALS, IF THERE ARE TWO
C                 ALPHAS OR TWO BETAS OF THE SAME IRREP
C                 KCRIT GUARDS AGAINS THE CASE OF ONE ALPHA AND ONE BETA
C                 WHEN THERE IS NO PHASE
                  NEQJMO=1
                  KSAVE=K
  505             CONTINUE
                  IF(IRRID(JMO).NE.IRRID(JMONEXT).OR.K.EQ.KCRIT)
     *                  GO TO 507
                     K=K+1
                     NEQJMO=NEQJMO+1
                     JMO=JMONEXT
                     IF(K.LT.NAEL) THEN
                        JMONEXT=ABS(IDTCSF(K+1,J))-NCORE
                        GO TO 505
                     END IF
  507             CONTINUE
C                 IF(IRRID(JMO).EQ.IRRID(JMONEXT).AND.K.NE.KCRIT) THEN
C                   FACT=FACT*(FF(IMOR,IT,JMO)*FF(IMORNEXT,IT,JMONEXT)-
C    *                         FF(IMOR,IT,JMONEXT)*FF(IMORNEXT,IT,JMO))
                  IF(NEQJMO.GT.1) THEN
                     FACT=FACT*DETFF(NEQJMO,FF,KSAVE,I,J,IT,NAEL,IDTCSF,
     *                               IRROW,NCORE)
C                   K=K+1
                    GO TO 510
                  END IF
                END IF
                FACT=FACT*FF(IMOR,IT,JMO)
  510           K=K+1
              IF(K.LE.NAEL) GO TO 500
              SUM1=SUM1+FACT
  600       CONTINUE
            SUM=SUM+SUM1*CDET(J)
C              WRITE(IW,*) 'SUMMING',J,SUM1,CDET(J)
  700       CONTINUE
            SUM=(SUM*LAMBD0(IG))/NT
            IF(ABS(SUM).GE.TINY.AND.DEBUG) WRITE(IW,9050) IRPNAM(IG),I,
     *                                                   SUM*SUM
            WEIGHT=WEIGHT+SUM*SUM
  800    CONTINUE
         CALL DDI_GSUMF(2309,WEIGHT,1)
         IF(WEIGHT.GT.SMALL) THEN
            IF(MASWRK) WRITE(IW,9060) IRPNAM(IG),WEIGHT
            NSYMG=NSYMG+1
            IRRSTAT=IOR(ISHFT(IRRSTAT,4),IG)
         END IF
         TOTWEIGHT=TOTWEIGHT+WEIGHT
         IF(NOIRR.LT.0.AND.TOTWEIGHT.GT.ONE-SMALL) THEN
            IF(DEBUG) WRITE(IW,*) 'SKIPPING THE REST OF IRRS',TOTWEIGHT
            GO TO 910
         END IF
  900 CONTINUE
C
  910 CONTINUE
      NBIT=64/NWDVAR
      IF(NSYMG.GT.NBIT/4) THEN
         WRITE(IW,*) 'SO MANY AN IRREP DO NOT FIT INTO ONE WORD!'
         IRRSTAT=0
      END IF
      RETURN
C
 9000 FORMAT(I7,F9.5,100I4)
 9050 FORMAT(1X,'THE PROJECTION ONTO ',A4,'CONTAINS A DETERMINANT ',I7,
     *       ' OF WEIGHT',1P,E11.4)
 9060 FORMAT(1X,'THE PROJECTION OF THIS CI STATE ONTO SPACE SYMMETRY ',
     *         A4,' WEIGHS',1P,E11.4/)
C9200 FORMAT(/1X,'POINT GROUPS WITH 3+ FOLD IRREPS ARE NOT READY YET.')
      END
C*MODULE MTHLIB  *DECK SEREP
      SUBROUTINE SEREP(NAO,NCORE,FF,CHARR,IRRID,IRROW,CLCAO,S,WORK1,
     *                 FUNSYM,PROJ,SMALL,PUREE,BAIL)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,DEBUG,PUREE,BAIL
C
      PARAMETER (MXSH=5000, MAXL=3, MXGRP=48, MXATM=2000,
     *           MXGTOT=20000, MXAO=8192)
      PARAMETER (ONE=1.0D+00,ZERO=0.0D+00)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SOSYM / EULANG(4,MXGRP),GAM(MXGRP,48),IRMON(MXAO)
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,NSAFMO
      COMMON /SYMREP/ IRPNAM(14),IPA(14),LAMBDA(14),LAMBD0(14),
     *                IADDR1(14),IADDR2(14),IADDR3(14)
      COMMON /SYMSPD/ PTR(3,3,48),DTR(6,6,48),
     *                FTR(10,10,48),GTR(15,15,48)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
      DIMENSION FF(MAXL,NT,NAO),CHARR(NT,NIRRED),IRRID(NAO),IRROW(NAO),
     *          CLCAO(NUM,NQMT),S(*),WORK1(NUM),FUNSYM(NUM,MAXL),
     *          PROJ(NUM,MAXL,MAXL),JOLT(MAXL),ALP(MAXL)
      DATA DBGME/8HIRRSTATE/,DEBUGG/8HDEBUG   /,SALC/8HSALC    /
C
C     THIS ROUTINE PREPARES SYMMETRY RELATED INFO TO BE USED FOR
C     NON-ABELIAN GROUPS (PUREE=.FALSE.)
C     PUREE=.TRUE. COOKS PUREE OUT OF SYMMETRY CONTAMINATED ORBITALS
C     IN THIS CASE A GUESS FROM SYMMOS IS USED ALTHOUGH IT IS POSSIBLE
C     TO GRIND THE ORBITALS AND MAKE PURE OUT OF THAT
C     NOTE THAT -NAO- STANDS FOR NUMBER OF ACTIVE ORBITALS
C     IF BAIL IS .TRUE. THEN THE CODE WOULD ABORT IF PROJECTION OF
C     AN ORBITAL ONTO SOME IRRED IS ZERO. BAIL=.FALSE. MAKES SENSE ONLY
C     FOR ORBITAL PURIFICATION WHEN SYMMOS GENERATED INCORRECT GUESS
C     (E.G., ASSIGNED B1,B1 FOR PX+PY,PX-PY ORBITALS IN C2V)
C     IF BAIL=.FALSE. AND THE ORBITALS DO NOT HAVE NON-ZERO PROJECTION
C     ONTO THE SPACE OF THE IRREDS SET BY IRMON, BAIL IS SET TO .TRUE.
C
      DEBUG=(EXETYP.EQ.DBGME.OR.EXETYP.EQ.DEBUGG.OR.EXETYP.EQ.SALC)
     *      .AND.MASWRK
C     IF(MASWRK.AND.LINEAR.AND.ICI.EQ.1) WRITE(IW,9220)
      INDEX=1
C     FIND CHARACTERS OF THE IRREPS
      DO 100 IG=1,NIRRED
         DO 200 IT=1,NT
            CHARR(IT,IG)=TRACEW(GAM(IT,INDEX),MXGRP,LAMBD0(IG))
  200    CONTINUE
         INDEX=INDEX+LAMBD0(IG)*LAMBD0(IG)
  100 CONTINUE
      I=1
      IID=0
C     ASSIGN A UNIQUE NUMBER TO EACH IRREP OF MOS.
  300 CONTINUE
         IID=IID+1
         IG=IRMON(I+NCORE)
         LAM=LAMBD0(IG)
         DO 400 J=1,LAM
            IF(IRMON(I+NCORE).NE.IG) THEN
               WRITE(IW,9300) I+NCORE,IRMON(I+NCORE),IG
               CALL ABRT
            END IF
            IRRID(I)=IID
            IRROW(I)=J
            I=I+1
  400    CONTINUE
      IF(I.LE.NAO) GO TO 300
C
      IF(DEBUG) WRITE(IW,9306) 'ACTIVE ORB IRRID',(IRRID(I),I=1,NAO)
      IF(DEBUG) WRITE(IW,9306) 'ACTIVE ORB IRROW',(IRROW(I),I=1,NAO)
 9306 FORMAT(1X,A,(1X,15I4/))
C
C     USE PROJECTION OPERATORS TO EXTRACT ROW-SYMMETRISED MOS
C
C     THEN FORM MATRICES OF "SEREDUCIBLE" REPRESENTATIONS FOR EACH AO
C     THIS MEANS TG X MO(I) = SUM OVER J OF FJ(G,I) SMO(J,I)
C     FJ(G,I) BEING SEREDUCIBLE REPRESENTATION
C     MO(I) ORIGINAL MO NUMBER I
C     TG GROUP SYMMETRY OPERATION
C     SMO ROW-SYMMETRISED MOS (ROW J FOR MO I)
C
C     LOOP OVER MOS
C
      DO 1000 I=1+NCORE,NAO+NCORE
        INDEX=1
        IG=IRMON(I)
        IGROW=IRROW(I-NCORE)
        LAM=LAMBD0(IG)
        DO 550 IG1=1,IG-1
          INDEX=INDEX+LAMBD0(IG1)*LAMBD0(IG1)
  550   CONTINUE
        IF(IGROW.EQ.1) CALL VICLR(JOLT,1,MAXL)
        DO 900 IL=1,LAM
          IAO=0
          IF(DEBUG) WRITE(IW,*) 'MO',I,' IRREP',IG,' ROW',IL
C         LOOP OVER AOS
          DO 800 II=1,NSHELL
            ITYP=KTYPE(II)
            IMIN=KMIN(II)
            IMAX=KMAX(II)
            KKEND=IMAX-IMIN
            DO 700 KK=0,KKEND
              SUM1=ZERO
              IAO=IAO+1
              DO 650 IT=1,NT
                KT=KLOC(MAPSHL(II,INVT(IT)))
                SUM=ZERO
                DO 600 KK1=0,KKEND
                  L=ITYP-1
                  M=KK+1
                  M1=KK1+1
                  IF(ITYP.EQ.2.AND.KKEND.EQ.3) THEN
C                   THE INFLAMOUS L-SHELLS
                    IF(KK.EQ.0.AND.KK1.NE.0.OR.KK.NE.0.AND.KK1.EQ.0)
     *                GO TO 600
                    L=1
                    M=KK
                    M1=KK1
                    IF(KK.EQ.0.AND.KK1.EQ.0) L=0
                  END IF
                  IF(L.EQ.0) F=ONE
                  IF(L.EQ.1) F=PTR(M,M1,IT)
                  IF(L.EQ.2) F=DTR(M,M1,IT)
                  IF(L.EQ.3) F=FTR(M,M1,IT)
                  IF(L.EQ.4) F=GTR(M,M1,IT)
                  SUM=SUM+CLCAO(KT+KK1,I)*F
600             CONTINUE
                SUM1=SUM1+SUM*GAM(IT,INDEX+(IL-1)*LAM+IL-1)
C               GAM IS TO BE COMPLEX-CONJUGATED
  650         CONTINUE
              SUM1=(SUM1*LAM)/NT
              IF(DEBUG) WRITE(IW,*) 'AO',IAO,' ORIG',CLCAO(KT+KK,I),
     *                               ' TRAN',SUM1
              PROJ(IAO,IGROW,IL)=SUM1
  700       CONTINUE
  800     CONTINUE
  900   CONTINUE
C
C       FIND ROW-SYMMETRISED NORMALISED FUNCTIONS.  IT IS IMPLIED
C       THAT ROWS ARE INCREASING SO THAT IGROW==LAM IS THE LAST
C
        IF(IGROW.EQ.LAM) THEN
          NJOLTED=0
          DO 930 IL=1,LAM
            DO 920 IL1=1,LAM
              CALL MTARBR(S,NUM,PROJ(1,IL,IL1),1,WORK1,NUM,1)
              ALPHA=DDOT(NUM,PROJ(1,IL,IL1),1,WORK1,1)
C
C               CHECK IF ORTHOGONAL TO ALL PREVIOUS
C             .09 CAN BE GENERALISED 1/MAXL/2 OR SO
C
              IF(ABS(ALPHA).GT.0.9D-01.AND.JOLT(IL1).EQ.0) THEN
                DO 910 MM=1,LAM
                  IF(JOLT(MM).EQ.0) GO TO 910
                  CALL MTARBR(S,NUM,PROJ(1,IL,IL1),1,WORK1,NUM,1)
                  ALPHA1=DDOT(NUM,FUNSYM(1,MM),1,WORK1,1)
                  IF(ABS(ALPHA1).GT.SMALL) GO TO 920
  910           CONTINUE
                NJOLTED=NJOLTED+1
                JOLT(IL1)=1
C               THOSE PSEUDO-ABELIAN GROUPS AMAZINGLY HAVE
C               IDENTICAL PROJECTORS FOR DIFFERENT ROWS (???)
                CALL DCOPY(NUM,PROJ(1,IL,IL1),1,FUNSYM(1,IL1),1)
                CALL DSCAL(NUM,ONE/SQRT(ALPHA),FUNSYM(1,IL1),1)
                IF(DEBUG) WRITE(IW,*) 'JOLT',IL1,' FROM ',IL,
     *                                (FUNSYM(JJ,IL1),JJ=1,NUM)
              END IF
  920       CONTINUE
  930     CONTINUE
C
          IF(NJOLTED.NE.LAM) THEN
            IF(BAIL) THEN
               WRITE(IW,9305) NJOLTED,LAM,I
               CALL ABRT
            ELSE
               BAIL=.TRUE.
               RETURN
            END IF
          END IF
C
          DO 960 ILL=1,LAM
            IMO=I-LAM+ILL-NCORE
            DO 935 IL=1,LAM
              CALL MTARBR(S,NUM,PROJ(1,ILL,IL),1,WORK1,NUM,1)
              ALP(IL)=DDOT(NUM,FUNSYM(1,IL),1,WORK1,1)
              IF(DEBUG) WRITE(IW,*) 'MULT',(PROJ(IUU,ILL,IL),IUU=1,NUM)
              IF(DEBUG) WRITE(IW,*) 'BY',(FUNSYM(IUU,IL),IUU=1,NUM)
C           ALPHA IS EXPANSION COEF OF MO IN TERMS OF SMO
  935       CONTINUE
            AA=DDOT(LAM,ALP,1,ALP,1)
            IF(DEBUG) WRITE(IW,*) 'AAAA',AA,(ALP(IUU),IUU=1,LAM)
            IF(MASWRK.AND..NOT.PUREE) WRITE(IW,9000)IMO+NCORE,IRPNAM(IG)
     *                                   ,(ALP(L)**2/AA*1.0D+02,L=1,LAM)
C           IF(ABS(AA-ONE).GT.TOL.AND.MASWRK) WRITE(IW,9100) AA
            IF(PUREE) THEN
               CALL DCOPY(NUM,FUNSYM(1,ILL),1,CLCAO(1,IMO+NCORE),1)
            ELSE
            DO 950 IT=1,NT
              DO 940 IL=1,LAM
                SUM=ZERO
                DO 937 IL1=1,LAM
                  SUM=SUM+GAM(IT,INDEX+(IL-1)*LAM+IL1-1)*ALP(IL1)
  937           CONTINUE
                FF(IL,IT,IMO)=SUM
  940         CONTINUE
             IF(DEBUG) WRITE(IW,*) 'IRRRS',I,IT,(FF(IL,IT,IMO),IL=1,LAM)
  950       CONTINUE
            END IF
  960     CONTINUE
        END IF
 1000 CONTINUE
      RETURN
C
 9000 FORMAT(1X,'ORBITAL',I3,'(',A4,')',' HAS ROW FRACTIONS ',
     *          5(F4.0,'% '))
 9300 FORMAT(1X,'THE ACTIVE ORBITALS MUST BE ORDERED IN SUCH WAY',
     *          ' THAT DEGENERATE'/
     *       1X,'MO''S ARE ADJACENT. THIS MO IS OUT OF ORDER',
     *           I3,' : ',2I2,/
     *       1X,'PLEASE REORDER THE ORBITALS. NOTE THAT',
     *          ' $GUESS IS NOT READ BY RUNTYP=TRANSITN.',/)
 9305 FORMAT(1X,'ERROR IN SEREP: NJOLTED,LAM=',2I9,' FOR MO=',I7/
     *       1X,'NOT ALL LINEARLY INDEPENDENT FUNCTIONS FOUND.')
      END
C*MODULE MTHLIB  *DECK TRACEW
      FUNCTION TRACEW(A,LDA,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(LDA,N*N)
C
C     RETURNS THE TRACE (SPUR) OF A SQUARE MATRIX NXN
C     W STANDS FOR WEIRD AS A MILD TERM FOR THE
C     WAY IRRED MATRICES ARE STORED
C
C     WARNING: NOT FOR GENERAL USE!!
C
      TRACEW=0.0D+00
      DO 100 I=1,N
         TRACEW=TRACEW+A(1,(N+1)*(I-1)+1)
  100 CONTINUE
      RETURN
      END
C*MODULE MTHLIB  *DECK IDETDISCO
      FUNCTION IDETDISCO(NCORE,NAEL,ICSF,JCSF,IRRID)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION ICSF(NAEL),JCSF(NAEL),IRRID(*)
C
      DO 100 K=1,NAEL
C
C        IN DETERMINING DISCOINCIDENCES IT IS ASSUMED THAT BOTH DETS
C        HAVE THE SAME MS, AS THEN IT IS GUARANTEED THAT BOTH HAVE
C        EQUAL NUMBER OF ALPHA'S AND BETA'S AND NO FALSE DISCOS DUE TO
C        SIGN OF ICSF (I.E. ALPHA-NESS OR BETA-NESS) ARISES
C
         IF(IRRID(ABS(ICSF(K))-NCORE).NE.IRRID(ABS(JCSF(K))-NCORE)) THEN
            IDETDISCO=0
            RETURN
         END IF
  100 CONTINUE
      IDETDISCO=1
      RETURN
      END
C*MODULE MTHLIB  *DECK LRIRREP
      SUBROUTINE LRIRREP(IRRL,IRRR)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      CHARACTER*1 LSYM(3)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGRP=48, MXAO=8192)
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00)
C
      DIMENSION IRRL(3),SUML(3),CHARL(3,MXGRP),
     *          IRRR(3),SUMR(3),CHARR(3,MXGRP)
C
      COMMON /SYMSPD/ PTR(3,3,48),DTR(6,6,48),
     *                FTR(10,10,48),GTR(15,15,48)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,NSAFMO
      COMMON /SYMREP/ IRPNAM(14),IPA(14),LAMBDA(14),LAMBD0(14),
     *                IADDR1(14),IADDR2(14),IADDR3(14)
      COMMON /SOSYM / EULANG(4,MXGRP),GAM(MXGRP,48),IRMON(MXAO)
      COMMON /SYMQMT/ IRPLAB(14),IRPNUM(14),IRPDIM(14),IRPDEG(14)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA LSYM/'X','Y','Z'/
C
C     THIS ROUTINE DETERMINES TO WHICH IRREPS THE THREE COMPONENTS
C     OF L (ANGULAR MOMENTUM) AND R (TRANSLATION) BELONG
C     PHASE COMES INTO PLAY BECAUSE L IS A PSEUDOVECTOR
C     (IE INVARIANT UNDER INVERSION WHEREAS R IS A VECTOR)
C     THE ORDER OF X,Y,Z IS AS IN PX,PY,PZ (1,2,3)
C     NOTE THAT LX CAN BELONG TO AT MOST TWO DIFFERENT REPRESENTATIONS
C     BUT ONLY IF COMPLEX BASIS SETS ARE USED, SO THAT IN GROUPS LIKE
C     C3 ONE HAS COMPLEX VALUED IRREPS.
C
      CALL VCLR(CHARL,1,3*NT)
      CALL VCLR(CHARR,1,3*NT)
      TOL=1.0D-10
      DO 100 IT=1,NT
         PHASE=ONE
         IF(DET3(T(1+(IT-1)*9)).LT.-TOL) PHASE=-ONE
         DO 100 I=1,3
            CHARL(I,IT)=CHARL(I,IT)+PHASE*PTR(I,I,IT)
            CHARR(I,IT)=CHARR(I,IT)+      PTR(I,I,IT)
  100    CONTINUE
      CALL VICLR(IRRL,1,3)
      CALL VICLR(IRRR,1,3)
      INDEX=1
      DO 200 IG=1,NIRRED
         CALL VCLR(SUML,1,3)
         CALL VCLR(SUMR,1,3)
         DO 300 IT=1,NT
            CHART=TRACEW(GAM(IT,INDEX),MXGRP,LAMBD0(IG))
            DO 400 I=1,3
               SUML(I)=SUML(I)+CHART*CHARL(I,IT)
               SUMR(I)=SUMR(I)+CHART*CHARR(I,IT)
C              COMPLEX CONJUGATE OF CHART
  400       CONTINUE
  300    CONTINUE
         DO 500 I=1,3
            SUML(I)=(SUML(I)*LAMBD0(IG))/NT
            SUMR(I)=(SUMR(I)*LAMBD0(IG))/NT
            IF(ABS(SUML(I)-ONE).GT.TOL.AND.ABS(SUML(I)).GT.TOL.OR.
     *         ABS(SUMR(I)-ONE).GT.TOL.AND.ABS(SUMR(I)).GT.TOL) THEN
            IF(ABS(SUML(I)-TWO).LT.TOL.AND.ABS(SUMR(I)-TWO).LT.TOL) THEN
C              MUST BE THE PSEUDO-ABELIAN GROUPS CN OR CNH OR DNH
            ELSE
              IF(MASWRK) WRITE(IW,*) 'UNEXPECTED SUMS. ',SUML(I),SUMR(I)
            END IF
            END IF
            IF(ABS(SUML(I)).GT.TOL) THEN
C              WRITE(IW,9000) 'L',LSYM(I),IRPLAB(IG),SUML(I)
               IF(IRRL(I).NE.0) THEN
                  IF(MASWRK) WRITE(IW,*) 'CANNOT ACCEPT TWO LABELS'
                  CALL ABRT
               END IF
               IRRL(I)=IG
            END IF
            IF(ABS(SUMR(I)).GT.TOL) THEN
C              WRITE(IW,9000) 'R',LSYM(I),IRPLAB(IG),SUMR(I)
               IF(IRRR(I).NE.0) THEN
                  IF(MASWRK) WRITE(IW,*) 'CANNOT ACCEPT TWO LABELS'
                  CALL ABRT
               END IF
               IRRR(I)=IG
            END IF
  500    CONTINUE
         INDEX=INDEX+LAMBD0(IG)*LAMBD0(IG)
  200 CONTINUE
      IF(MASWRK) WRITE(IW,8900)
      DO 600 I=1,3
         IF(IRRL(I).EQ.0) THEN
            IF(MASWRK) WRITE(IW,*) 'CANNOT IDENTIFY L',LSYM(I)
            CALL ABRT
         END IF
         IF(MASWRK) WRITE(IW,9000) 'L',LSYM(I),IRPLAB(IRRL(I))
  600 CONTINUE
      DO 700 I=1,3
         IF(IRRR(I).EQ.0) THEN
            IF(MASWRK) WRITE(IW,*) 'CANNOT IDENTIFY R',LSYM(I)
            CALL ABRT
         END IF
         IF(MASWRK) WRITE(IW,9000) 'R',LSYM(I),IRPLAB(IRRR(I))
  700 CONTINUE
      RETURN
 8900 FORMAT(1X,'THE SYMMETRIES OF THE ANGULAR MOMENTUM L AND ',
     *          'TRANSLATION R:')
 9000 FORMAT(1X,2A1,' BELONGS TO ',A4)
      END
C*MODULE MTHLIB  *DECK DETFF
      DOUBLE PRECISION FUNCTION DETFF(N,FF,K,I,J,IT,NAEL,
     *                               IDTCSF,IRROW,NCORE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000,MXSH=5000,MAXL=3)
      DIMENSION IDTCSF(NAEL,*),FF(MAXL,NT,*),IRROW(*),A(MAXL,MAXL)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
C     THIS IS A SPECIAL KIND OF MATRIX DETERMINANT
C     (MORE LIKE MINOR THAN DETERMINANT)
C
C                 IF(IRRID(JMO).EQ.IRRID(JMONEXT).AND.K.NE.KCRIT) THEN
C                   FACT=FACT*(FF(IMOR,IT,JMO)*FF(IMORNEXT,IT,JMONEXT)-
C    *                         FF(IMOR,IT,JMONEXT)*FF(IMORNEXT,IT,JMO))
C     WRITE(IW,*) 'ENTERING DETFF',N,I,J,K
      DO 100 II=1,N
         IMOR=IRROW(ABS(IDTCSF(K+II-1,I))-NCORE)
         DO 100 JJ=1,N
            JMO=ABS(IDTCSF(K+JJ-1,J))-NCORE
            A(II,JJ)=FF(IMOR,IT,JMO)
  100 CONTINUE
C     NOW ALL IS LEFT IS TO GET THE DETERMINANT
C     DETFF=DET(A)
      IF(N.EQ.2) DETFF=A(1,1)*A(2,2)-A(1,2)*A(2,1)
      IF(N.EQ.3) DETFF=DET3(A)
C     DET3 REQUIRES LEADING DIMENSION 3
      IF(N.NE.2.AND.N.NE.3) THEN
         WRITE(IW,*) 'DET FAILED'
         CALL ABRT
      END IF
C     WRITE(IW,*) 'EXITING DETFF'
      RETURN
      END
C*MODULE MTHLIB  *DECK CSFTODET
      SUBROUTINE CSFTODET(NORB,CCSF,MXSPIN,NAEL,MAXDET,NDET,IDTCSF,
     *                    CDET,ICASE,IICAS,IPRIM,CGC,BUFFPK,NTRAP,BK)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXAO=8192, MXATM=2000)
      PARAMETER (ONE=1.0D+00)
      DIMENSION IDTCSF(NAEL,*),CDET(MAXDET),BUFFPK(MXSPIN,*),ICASE(*),
     *          IPRIM(MXSPIN,*),CGC(*),IORDER(MXAO),IICAS(MXSPIN),
     *          BK(MXSPIN),NTRAP(MXSPIN),ICSF(MXAO),N1234(4)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /ORBSET/ NORBMX,NORBS,NCORBS,NLEVS,NA,NB,NC,NSYM,MSYM,
     *                IDOCC,IVAL,IMCC,ISYM(MXAO),ICODE(MXAO),
     *                NLCS(MXAO),LEVPT(MXAO),LEVNR(MXAO),IOUT(MXAO),
     *                NREFS,IEXCT,NFOCI,INTACT
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,MA,MB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      DATA N1234/ 0, 1, -1, 0 /
C
C     DATA DEBUG/8HDEBUG   /, DBUGME/8HCSF2DET /
C
C     NORBS
      DO 770 I=1,NORB
         DO 765 J=1,NORBMX
            IF(IOUT(J).EQ.I) IORDER(I)=J
  765    CONTINUE
  770 CONTINUE
C
      IPX = 0
      IBK = 0
      NPHI = 1
C     NORBS
      DO 40 J=1,NORB
C           WRITE(IW,*) 'NOOO',J,ICASE(J)
         IBK = IBK + N1234(ICASE(J))
         GO TO (40,42,42,44), ICASE(J)
C
   42    CONTINUE
         IPX=IPX+1
         IICAS(IPX)=ICASE(J)
         GO TO 40
C
   44    CONTINUE
         ICHK = IBK/2
         ICHK = IBK - 2*ICHK
         IF(ICHK.EQ.1) NPHI =-NPHI
   40 CONTINUE
C
      NSORB = 0
      CALL SPNFNC(IPX,MUL,IICAS,MXSPIN,IPRIM,CGC,NPRI,BUFFPK,NTRAP,BK)
C
      IF(NPHI.LT.0) CALL DSCAL(NPRI,-ONE,CGC,1)
      IF(NPRI.GE.2**MXSPIN) THEN
         IF(MASWRK) WRITE(IW,9020)
         CALL ABRT
      END IF
      DO 200 IPR=1,NPRI
        CALL MAKDET(NORB,ICASE,IPRIM(1,IPR),IORDER,ICSF,NSORB)
C       CALL MAKDET(NORBS,ICASE,IPRIM(1,IPR),IORDER,ICSF,NSORB)
C       THE CSFS HAVE NOT BEEN SORTED YET
        CALL DETREOR(NSORB,ICSF,NPHASE)
        IF(NPHASE.LT.0) CGC(IPR)=-CGC(IPR)
        CALL IFINDET(IRET,ICSF,IDTCSF,MAXDET,NDET,NAEL)
C          WRITE(IW,*) 'UPDATING',IRET,IPR,MAXDET,KSTATE,
C    *                 CGC(IPR)*CCSF(IWKS,KSTATE)
          CDET(IRET)=CDET(IRET)+CGC(IPR)*CCSF
  200 CONTINUE
      RETURN
 9020 FORMAT('THE NUMBER OF SLATER DETS WILL NOT FIT INTO 12 BITS')
      END
C*MODULE MTHLIB  *DECK KOMVEC
      FUNCTION KOMVEC(N,EPS,DA,DB,INCX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION DA(INCX,N),DB(INCX,N)
C
C     COMPARES TWO VECTORS, RETURNS THE INDEX OF THE FIRST DIFFERENT
C     ELEMENT, OR 0 IF EQUAL
C
      DO 50 I = 1,N
         IF(ABS(DA(1,I)-DB(1,I)).GT.EPS) THEN
            KOMVEC=I
            RETURN
         END IF
   50 CONTINUE
      KOMVEC=0
      RETURN
      END
C*MODULE MTHLIB  *DECK POLFIT
      SUBROUTINE POLFIT(NPTS,NTERMS,MODE,X,Y,A,SIGMAY,CHISQR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAXPOLY=10)
      DIMENSION X(NPTS),Y(NPTS),A(NTERMS),SIGMAY(NPTS)
      DIMENSION SUMX(19),SUMY(MAXPOLY),ARRAY(MAXPOLY,MAXPOLY)
      DATA ZERO/0.0D+00/ , ONE/1.0D+00/ , TWO/2.0D+00/
C
C     LEAST SQUARES FIT TO DATA WITH A POLYNOMIAL CURVE.
C     Y=A(1)+A(2)*X+A(3)*X**2+A(4)*X**3+...
C
C     X        DATA PTS FOR INDEPEPENDENT VARIABLE
C     Y        DATA PTS FOR DEPENDENT VARIABLE
C     SIGMAY   STD DEV FOR Y DATA PTS. (ONLY IF -MODE- IS 1)
C     NPTS     # OF PAIRS OF DATA POINTS X,Y
C     NTERMS   # OF COEFFICIENTS (DEGREE OF POLY+1). (<= MAXPOLY)
C     MODE      1   WEIGHT = 1/SIGMAY
C               0   WEIGHT = 1.0
C              -1   WEIGHT = 1/ABS(Y)
C     A        ARRAY OF COEFFICIENTS OF POLYNOMIAL (OUTPUT)
C     CHISQR   REDUCED CHI SQUARE FOR FIT'S VARIANCE. (OUTPUT)
C
C           POLFIT WAS ADAPTED FROM THE VERSION GIVEN IN
C     DATA REDUCTION AND ERROR ANALYSIS FOR THE PHYSICAL SCIENCES
C                    PHILIP R. BEVINGTON
C          MCGRAW-HILL BOOK CO., N.Y., 1969., PAGES 140-142.
C
      IF(NPTS.GT.20  .OR.  NTERMS.GT.10) THEN
         WRITE(6,*) 'POLFIT: BOGUS POLYNOMIAL FIT ARGUMENTS'
         CALL ABRT
      END IF
C
C     ACCUMULATE WEIGHTED SUMS
C
      NMAX=2*NTERMS-1
      DO 13 N=1,NMAX
   13 SUMX(N)=ZERO
      DO 15 J=1,NTERMS
   15 SUMY(J)=ZERO
      CHISQ=ZERO
C
      DO 50 I=1,NPTS
      XI=X(I)
      YI=Y(I)
C
      IF(MODE) 32,37,39
   32 IF(YI) 35,37,33
   33 WEIGHT=ONE/YI
      GO TO 41
   35 WEIGHT=ONE/(-YI)
      GO TO 41
   37 WEIGHT=ONE
      GO TO 41
   39 WEIGHT=ONE/(SIGMAY(I)*SIGMAY(I))
C
   41 XTERM=WEIGHT
      DO 44 N=1,NMAX
      SUMX(N)=SUMX(N)+XTERM
   44 XTERM=XTERM*XI
      YTERM=WEIGHT*YI
      DO 48 N=1,NTERMS
      SUMY(N)=SUMY(N)+YTERM
   48 YTERM=YTERM*XI
      CHISQ=CHISQ+WEIGHT*YI*YI
   50 CONTINUE
C
C     CONSTRUCT MATRICES AND CALCULATE COEFFICIENTS.
C
      DO 54 J=1,NTERMS
      DO 54 K=1,NTERMS
      N=J+K-1
   54 ARRAY(J,K)=SUMX(N)
      CALL POLDET(DELTA,ARRAY,NTERMS,MAXPOLY)
C
      IF(DELTA) 61,57,61
   57 CHISQR=ZERO
      DO 59 J=1,NTERMS
   59 A(J)=ZERO
      GO TO 80
C
   61 DO 70 L=1,NTERMS
      DO 66 J=1,NTERMS
      DO 65 K=1,NTERMS
      N=J+K-1
   65 ARRAY(J,K)=SUMX(N)
   66 ARRAY(J,L)=SUMY(J)
      CALL POLDET(DET,ARRAY,NTERMS,MAXPOLY)
      A(L)=DET/DELTA
   70 CONTINUE
C
C     CALCULATE CHI SQUARE
C
      DO 75 J=1,NTERMS
      CHISQ=CHISQ-TWO*A(J)*SUMY(J)
      DO 75 K=1,NTERMS
      N=J+K-1
   75 CHISQ=CHISQ+A(J)*A(K)*SUMX(N)
      FREE=NPTS-NTERMS
      IF(FREE.EQ.ZERO) FREE=ONE
      CHISQR=CHISQ/FREE
   80 RETURN
      END
C*MODULE MTHLIB  *DECK POLDET
      SUBROUTINE POLDET(DET, ARRAY, NORDER, MAXPOLY)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION  ARRAY (MAXPOLY,MAXPOLY)
      DATA ZERO/0.0D+00/, ONE/1.0D+00/
C
C     THIS FUNCTION CALCULATES THE DETERMINATE OF A SQUARE MATRIX.
C
C     ARRAY    THE MATRIX OF ORDER NORDER WHICH IS TO BE EVALUATED.
C              THIS SUBPROGRAM DESTROYS THE MATRIX ARRAY
C     NORDER   THE ORDER OF THE SQUARE MATRIX TO BE EVALUATED.
C
C     FUNCTION SUBPROGRAM DETERM WAS ADAPTED FROM THE SAME SOURCE AS
C     POLFIT. DETERM APPEARED IN THE FIRST APPENDIX OF THAT BOOK.
C
      DET=ONE
      DO 50 K=1,NORDER
C
C     INTERCHANGE COLUMNS IF DIAGONAL IS ZERO
C
      IF (ARRAY(K,K))41,21,41
 21   DO 23 J=K,NORDER
      IF(ARRAY(K,J))31,23,31
 23   CONTINUE
      DET=ZERO
      GO TO 60
  31  DO 34 I=K,NORDER
      SAVE=ARRAY(I,J)
      ARRAY(I,J)=ARRAY(I,K)
  34  ARRAY(I,K)=SAVE
      DET=-DET
C
C     SUBTRACT ROW K FROM LOWER ROWS TO GET DIAGONAL MATRIX
C
 41   DET=DET*ARRAY(K,K)
      IF(K-NORDER)43,50,50
 43   K1=K+1
      DO 46 I=K1,NORDER
      DO 46 J=K1,NORDER
 46   ARRAY(I,J)=ARRAY(I,J)-ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
 50   CONTINUE
 60   RETURN
      END
C*MODULE MTHLIB  *DECK VECROT
      SUBROUTINE VECROT(V1,V2,A)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION V1(3),V2(3),A(3,3),V12(3)
      PARAMETER(ONE=1.0D+00,SMALL=1.0D-12)
C
C     GENERATE A 3X3 UNITARY MATRIX THAT ROTATES VECTOR V1 TO V2.
C     THE VECTORS NEED NOT BE NORMALISED.
C
      VN1= V1(1)*V1(1)+V1(2)*V1(2)+V1(3)*V1(3)
      VN2= V2(1)*V2(1)+V2(2)*V2(2)+V2(3)*V2(3)
C     IF EITHER IS A NULL VECTOR SET THE MATRIX TO BE UNIT
      IF(VN1.LT.SMALL.OR.VN2.LT.SMALL) THEN
         CALL RUNITV(3,3,A)
         RETURN
      ENDIF
C
C     THE ORDER V1 AND V2 BELOW IS SUCH THAT THE DEFINITION OF THE
C     VECTOR PRODUCT AND THE ROTATION DIRECTION (ANTICLOCKWISE) AGREE.
C
      CALL VECPRD(V12,V2,V1)
      V12S=V1(1)*V2(1)+V1(2)*V2(2)+V1(3)*V2(3)
      COSW=V12S/SQRT(VN1*VN2)
C     IF V1 IS COLLINEAR TO V2: UNIT MATRIX
      IF(ABS(COSW-ONE).LT.SMALL) THEN
         CALL RUNITV(3,3,A)
         RETURN
      ENDIF
      IF(ABS(COSW+ONE).LT.SMALL) THEN
C
C     IF V1 IS ANTICOLLINEAR TO V2, V12 IS A (NEARLY) ZERO VECTOR.
C       IN THIS CASE THERE IS AN INIFINITE NUMBER OF VV12 (LYING IN
C       THE SAME PLANE ORTHOGONAL TO BOTH V1 AND V2). A HARD-WIRED
C       CHOICE IS MADE THEN BY ORTHOGONALISING Z- OR Y-AXIS ORTH,
C       WHICHEVER IS LINEARLY INDEPENDENT.
C
        IF(ABS(V1(3)).GT.SMALL) THEN
          V12(1)=VN1*V1(2)*V1(1)
          V12(2)=VN1*V1(2)*V1(2)-ONE
          V12(3)=VN1*V1(2)*V1(3)
        ELSE
          V12(1)=VN1*V1(3)*V1(1)
          V12(2)=VN1*V1(3)*V1(2)
          V12(3)=VN1*V1(3)*V1(3)-ONE
        ENDIF
      ENDIF
      VV12=SQRT(V12(1)*V12(1)+V12(2)*V12(2)+V12(3)*V12(3))
      V12(1)=V12(1)/VV12
      V12(2)=V12(2)/VV12
      V12(3)=V12(3)/VV12
      OMEGA=ACOS(COSW)
      CALL MATNOM(OMEGA,V12,A)
      RETURN
      END
C*MODULE MTHLIB  *DECK MATNOM
      SUBROUTINE MATNOM(OMEGA,N,A)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION N
      PARAMETER(ONE=1.0D+00)
      DIMENSION N(3),A(3,3)
C
C     GENERATE A 3X3 ROTATION MATRIX FROM THE N,OMEGA REPRESENTATION
C     WHERE N IS AN ORTH DEFINING A ROTATION AXIS AND OMEGA IS THE
C     ROTATION ANGLE. N IS AN ORTH AND AS SUCH MUST BE NORMALISED.
C     THIS ROUTINE GENERATES PURE ROTATION (EVEN PARITY).
C
      CO=COS(OMEGA)
      CO1=ONE-CO
      SI=SIN(OMEGA)
      A(1,1)=CO+CO1*N(1)*N(1)
      A(1,2)=CO1*N(1)*N(2)-SI*N(3)
      A(1,3)=CO1*N(1)*N(3)+SI*N(2)
      A(2,1)=CO1*N(2)*N(1)+SI*N(3)
      A(2,2)=CO+CO1*N(2)*N(2)
      A(2,3)=CO1*N(2)*N(3)-SI*N(1)
      A(3,1)=CO1*N(3)*N(1)-SI*N(2)
      A(3,2)=CO1*N(3)*N(2)+SI*N(1)
      A(3,3)=CO+CO1*N(3)*N(3)
      RETURN
      END
