C 24 MAR 07 - MWS - PAD FRGINF COMMON BLOCK
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C  7 SEP 06 - DMC - GVBITR: IMPLEMENT FULL SVPE SOLVATION
C 13 MAR 06 - MWS - ASSIGN GVB LAGRANGIAN TO UNIQUE DAF RECORD
C 14 NOV 05 - DGF,HL - SOLVENT SCF MEMORY CHANGES
C 19 SEP 05 - MWS - ADJUSTMENT FOR H OR I FUNCTION RUNS
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 27 JUN 05 - MWS - GVBITR: SVP MEMORY ALLOCATION CHANGE
C 30 APR 05 - MWS - IMPLEMENT SVP SOLVATION MODEL
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL, FRGINF
C  7 SEP 04 - MWS - PAD COMMON BLOCK INTFIL
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C  7 APR 04 - HL  - IEF-PCM/EFP FOR GVB
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  3 SEP 03 - MWS - GVBITR: INCLUDE MOPAC PEPTIDE BOND ENERGY
C 12 DEC 02 - MWS - MPCJK:NO DAREAD IF NO 2E- INTS, HSTARG:LOWER CUTOFF
C 22 MAY 02 - MWS - FIX PARALLEL ACAVO RUNS
C 17 APR 02 - KRG - GVBITR: ADD DUMMY ARG TO CALL FOR DIIS
C 24 JAN 02 - MWS - GENERATE ACAVO ORBITALS
C 30 NOV 01 - MWS - GVBITR: WIDEN ITERATION ENERGY FIELD
C 25 JUN 01 - MWS - ALTER COMMON BLOCK SCFOPT
C 28 JUL 00 - MWS - GVBCVG: SIMPLIFY TO REMOVE IBM OPTIMIZER CONFUSION
C 11 JUN 00 - MWS - GVBRHF: ALWAYS ALLOCATE A HONDO INTEGRAL BUFFER
C 16 FEB 00 - VK  - GVBFCK: ENABLED MO FREEZING IN GVB
C 13 MAR 99 - MWS - GVBITR: CORRECT CALL TO CLENMO ROUTINE
C 13 FEB 99 - MWS - FIX CALLS TO PUSQL AND TO MVOS
C 12 NOV 98 - GDF - CHANGE BIT PACKING TO ISHIFT
C 27 OCT 98 - MAF - GVBLAG,GVBCVG,GVBDIA,GVBFCK,GVBITR,GVBOUT,GVBTRN,
C                   RHFGVB: ALLOW FOR USE OF SPHERICAL HARMONICS
C 13 APR 98 - MWS - GVBOUT: SAVE SYMMETRY LABELS ON DAF
C  6 JAN 98 - MWS - CHANGE CALL TO RACLOS
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C 24 JUN 97 - MWS - GVBITR: CLOSE DIIS RANDOM FILE BEFORE EXIT
C 13 MAY 97 - MWS - GVBOUT: ALWAYS PRINT ASYMMETRY OF THE LAGRANGIAN
C 14 FEB 97 - MWS - GVBFCK: CHANGE VIRTUAL BLOCK WHEN NO CORES PRESENT
C 17 JAN 97 - MWS - GVBITR: GENERATE MVO-S AFTER SETTING UP LAGRANGIAN
C 18 DEC 96 - MWS - ADJUST CALL TO TWOEI
C 13 NOV 96 - DGF - ADD SEMIEMPIRICAL GVB
C 12 NOV 96 - MWS - GVBITR: ALLOW DIRECT MVO GENERATION
C 13 JUN 96 - VAG - ADD VARIABLE FOR CI TYPE TO SCFOPT COMMON
C 24 MAY 96 - GMC - GVBITR: ONLY SET NSO TO ZERO FOR PP(1) CASE.
C 11 APR 96 - GMC - CHANGE CALL TO SONEWT, AND WHEN SOSCF IS STARTED
C 14 SEP 95 - GMC - GVBHES: CHANGE DIAGONAL HESSIAN TOLERANCE
C 26 JUL 95 - GMC - IMPLEMENT SOSCF FOR ALL CASES, AND DIIS FOR 0/1 PAIR
C 21 APR 95 - MWS - USE CORRECT NAME FOR EXCHANGE INTEGRAL SCREENING
C  6 APR 95 - MWS - GVBITR: REMOVE AMPERSAND CONTINUATIONS
C  5 MAR 95 - MWS - GVBITR: CONVERGENCE IN OPT RUNS DEPENDS ON GRADIENT
C  1 FEB 95 - WC  - GVBITR: CHANGE ARGS TO EXTRAP,MVOS FOR MOROKUMA
C 12 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C  5 AUG 94 - MWS - HSTARG: DOUBLE LABEL PACKING, CHANGE TWOEI ARGS
C  1 AUG 94 - SPW - RHFGVB,GVBITR: ADD ABFLD TO CALL TO EFPCMP
C  1 JUN 94 - MWS - CHANGE JKUNIT NUMBER
C 27 MAY 94 - SPW - RHFGVB,GVBITR: CALC. IND. DIPOLE ON FRAGS EACH ITER
C 19 APR 94 - MWS - GVBFCK,GVBHAM,GVBITR,GVBLAG,HSTARG,NEWCI,OPTCI,ROTB:
C                   CHANGES TO PERMIT CALCULATION WITH NCO=0
C  1 APR 94 - SPW - GVBITR: INDUCE SCRF EVERY ITERATION
C  8 JUN 92 - TLW - PARALLELIZE DISK METHOD
C 24 MAR 92 - MWS - RENAME SELECT TO SELGV1 TO AVOID TCGMSG CONFLICT
C 21 MAR 92 - MWS - ALLOCATE ADDITIONAL MEMORY FOR DIRECT GVB
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C 12 MAR 92 - TLW - PARALLELIZE DIRECT METHOD
C  5 MAR 92 - MWS - GVBOUT: USE SYMMOS, PUT OLD LEVEL SHIFTER HERE
C 11 JAN 92 - TLW - MAKE PARALLEL I/O CHANGES
C 17 DEC 91 - MWS - USE FUZZIER CONV. CRIT. FOR DIRECT SCF RUNS
C 31 OCT 91 - JHJ - GVBITR: REMOVE ENUC CALL.
C 11 OCT 91 - MWS - GVBFCK,GVBITR,HSTARG,JKFORM,RHFGVB: CHANGES
C                   FOR DIRECT GVB, AND SIMPLIFIED DIIS ROUTINE.
C  7 OCT 91 - MWS - CHANGE UNIX BYTE UNPACKING
C 19 AUG 91 - JAM - GVBFCK: TRANSFORM ERROR MATRIX TO AO BASIS
C  7 JUL 91 - MWS - GVBITR: CHANGE USE OF ORTHO,QMATRX
C 29 JUN 91 - MWS - JKFORM: CALL SYMH FOR CORE,OPEN,PAIR OPERATORS
C 29 APR 91 - MWS - HSTARG: USE FULL 16 BIT UNPACKING ON PK FILES
C 30 OCT 90 - MWS - USE CLENMO ONLY FOR ABELIAN GROUPS
C 12 SEP 90 - MWS - INTRODUCE MXATM
C  6 JUL 90 - MWS - GVBTCH,NEWCI: SAVE PART OF CI HAMLTN IN /GVBCI/
C  3 JUN 90 - JAM - CHANGE VAX BYTE PACKING TO FULL 8 BITS
C  1 MAY 90 - MWS - FLUSH PRINT BUFFER ON EVERY ITERATION
C 21 DEC 89 - MWS - MORE DIGITS IN ENERGY WHEN PUNCHING MO-S
C 16 DEC 89 - MWS - HSTARG: TYPO FIXED, "ALL" SHOULD HAVE BEEN NJK.
C  3 DEC 89 - MWS - CHANGE FINAL PRINTOUT TO SHOW OCCUPATION NUMBERS,
C                   REMOVE A TWO FROM PAIR ORBITAL'S ENERGIES.
C  9 MAY 89 - MWS - REDDEN: GENERATE CORRECT A+B DM1 (MOST OF THE TIME)
C 20 FEB 89 - STE - HSTARG: FIX PREAD CALL (NXX WAS TYPED NXP)
C  1 FEB 89 - MWS - USE SYMTRZ, PREVS IN GVBOUT
C 15 NOV 88 - KF  - DIMENSION IWRK(L1) IN GVBITER, H2(L2) IN VIRSET
C 14 NOV 88 - MWS - GENERATE MVO-S AFTER GVB CONVERGENCE
C 24 AUG 88 - MWS - CORRECT NCO.EQ.0 CALL TO DCOPY IN JKFORM
C  9 JUL 88 - MWS - NEW GVBFCK FROM JOHN, AVOID CONVENTIONAL METHODS
C                   OF CONVERGENCE WHEN DIIS COMES ON.
C  3 JUL 88 - JAM - USE DIIS FOR THE GVB-PP(1) CASE.  ADD GVBFCK
C 30 MAY 88 - MWS - UP AO-S FROM 255 TO 2047
C  9 MAY 88 - MWS - READ ORBS FROM DAF BEFORE CALLING GVBOUT
C  7 MAY 88 - MWS - RENAME HSGFM TO JKMATS, BETTER COMMENTS
C 24 APR 88 - MWS - MASSIVE CLEANUP: MEMORY ALLOCATION DOES NOT
C                   USE SLOW MEMORY, USE EXETYP=GVB FOR DEBUG PRINT,
C                   ROUTINES DENGVL, DENGVB, GVBDAG, GVBDRV DELETED,
C                   SCFWFN AND GVBORB COMBINED IN COMMON GVBWFN
C 10 APR 88 - MWS - DON'T TAMPER WITH CLOSED OR OPEN SHELL EIGENVALUES
C 31 MAR 88 - MWS - PRINT TOTAL SCF TIME, CLEAN UP CONVERGENCE LOGIC,
C                   USE CLENMO DURING THE ITERATIONS.
C 13 NOV 87 - STE - USE EXETYP
C 28 OCT 87 - STE - GVBOUT: CHANGE 8010 TO IN LINE FORMAT
C 16 OCT 87 - MWS - PUNCH ENERGY AND TITLES WITH ORBITALS
C  6 AUG 87 - MWS - INCLUDE ETA VERSION
C 21 MAY 87 - MWS - ICH FROM 100 TO 255 IN SELGV1, BIT RSHIFT USED
C 28 APR 87 - STE - CHANGE CALLS TO GLDIAG(IN NEWCI), QMATRX, GVBDAG
C 24 APR 87 - STE - REDDEN: REMOVE IA
C                   GVBITR: REMOVE IA FROM CALLS TO REDDEN,ORTHO
C                   JKFORM: REMOVE IA FROM CALLS TO TFTRI
C 24 NOV 86 - MWS - CHANGE RSHIFT TO RRSHFT
C  3 NOV 86 - STE - SCHORT:USE DSCAL; UPOCBS:MOVE FAC=0; PARAMETERS
C                   REMOVE NAV REFERENCES
C  8 JUL 86 - JAB - ADD BIT PACKING/UNPACKING FOR CRAY, CELERITY
C  7 JUL 86 - JAB - RENAME ROUTINE "SHIFTR" TO "LVLGVB"
C  7 JUL 86 - MWS - SANITIZE FLOATING POINT CONSTANTS
C  6 NOV 85 - STE - DELETE /GVBFIL/, PUT JKUNIT IN /IOFILE/ AS IPK=9
C  5 NOV 85 - STE - USE GENERIC ABS,MAX,SQRT; CONVERGE WITH DAMP ON
C 27 JUL 85 - MWS - CHANGE IBM INDEX PACKING
C 23 APR 85 - MWS - FIX TYPO IN GVBOUT, REMOVE NEWORD ARG FOR GVBOUT
C 16 APR 85 - MWS - ONLY SORT VIRTUALS, FIX OPTCI RECORD LENGTH (DFF),
C                   PRINT ITERS WITH FINAL E
C 10 MAR 85 - MWS - CALL OPENJK/CLOSJK TO OPEN/CLOSE JKFILE
C  5 JAN 83 - STE - DELETE REWIND IPK'S, MOVE GVB BANNER TO GVBDRV
C 27 DEC 83 - STE - ORDER ORBITALS IN GVBOUT BEFORE DOING OUTPUT
C                   USE SQREAD/SQWRIT INSTEAD OF IMPLICIT LOOPS
C 26 DEC 83 - STE - MAKE GVBITR PRINT FINAL E AND CI COEFS ALWAYS
C 14 DEC 83 - STE - RENAME COMMON BLOCK DMPTST TO DMPING
C 18 NOV 83 - STE - FPS UNPACKING IN HSTARG
C  3 NOV 83 - STE - REMOVE BRANCH TO 160 IN RHFGVB TO DEFINE SM VARS
C  4 OCT 82 - MWS - CONVERT FOR IBM
C
C*MODULE GVB     *DECK CIEXPR
      SUBROUTINE CIEXPR(KCORB,IA,L1)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      DIMENSION IA(L1),KCORB(2,12)
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00, TM6=1.0D-06)
C
C     ----- UPDATE F, ALPHA, AND BETA MATRICES BASED ON
C           THE NEW CI COEFFICIENTS -----
C
      TEST = TM6
C
C     ----- NORMALIZE CI COEFFICIENTS AND UPDATE F -----
C
      DO 140 KPAIR = 1,NPAIR
         Y = ZERO
         DO 100 KORB = 1,2
            Y = Y+CICOEF(KORB,KPAIR)**2
  100    CONTINUE
         IF (Y .LT. TEST) Y = TEST
         Y = SQRT(Y)
         DO 120 KORB = 1,2
            CICOEF(KORB,KPAIR) = CICOEF(KORB,KPAIR)/Y
            KK = KCORB(KORB,KPAIR)
            NX = NCONF(KK)
            F(NX) = CICOEF(KORB,KPAIR)**2
  120    CONTINUE
  140 CONTINUE
C
C     ----- UPDATE ALPHA AND BETA -----
C
      DO 180 KPAIR = 1,NPAIR
         KONE = KCORB(1,KPAIR)
         KTWO = KCORB(2,KPAIR)
         LO = NCONF(KONE)
         LT = NCONF(KTWO)
         LOLO = IA(LO)+LO
         LOLT = IA(LT)+LO
         LTLT = IA(LT)+LT
         DO 160 K = 1,NHAM
            LOK = IA(LO)+K
            IF (K .GT. LO) LOK = IA(K)+LO
            LTK = IA(LT)+K
            IF (K .GT. LT) LTK = IA(K)+LT
            ALPHA(LOK) = TWO*F(LO)*F(K)
            ALPHA(LTK) = TWO*F(LT)*F(K)
            BETA(LOK) = -F(LO)*F(K)
            BETA(LTK) = -F(LT)*F(K)
  160    CONTINUE
         ANN = CICOEF(1,KPAIR)*CICOEF(2,KPAIR)
         ALPHA(LOLO) = F(LO)
         ALPHA(LTLT) = F(LT)
         BETA(LOLO) = ZERO
         BETA(LTLT) = ZERO
         BETA(LOLT) = ANN
         ALPHA(LOLT) = ZERO
  180 CONTINUE
      RETURN
      END
C*MODULE GVB     *DECK DENSHL
      SUBROUTINE DENSHL(TRANS,D,IHAM,L1,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      DIMENSION TRANS(L1,*),D(L2)
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
      PARAMETER (ZERO=0.0D+00)
C
C      ----- FORM DENSITY FOR THE IHAM'TH SHELL -----
C      ONLY ORBITALS OCCUPIED IN THIS SHELL MAKE A CONTRIBUTION
C
      IJ = 0
      DO 160 I = 1,L1
         DO 140 J = 1,I
            IJ = IJ+1
            SUM = ZERO
            DO 120 KX = 1,NORB
               IF(NCONF(KX).EQ.IHAM) SUM = SUM+TRANS(I,KX)*TRANS(J,KX)
  120       CONTINUE
            D(IJ) = SUM
  140    CONTINUE
  160 CONTINUE
      RETURN
      END
C*MODULE GVB     *DECK GVBCVG
      SUBROUTINE GVBCVG(HOLD,E,L1,L0,SQCDF,OUT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL OUT
C
      DIMENSION HOLD(L1,*),E(*)
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- CALCULATE THE SUM OF THE SQUARES OF THE OFF
C           DIAGONAL COEFFICIENTS IN THE TRANSFORMED BASIS
C           EXCLUDING TERMS CONNECTING MO'S IN THE SAME
C           FOCK OPERATOR -----
C
      SQCDF = ZERO
      DO 140 J = 1,NORB
         NCJ = NCONF(J)
         DO 100 I = 1,L0
            IF(NCONF(I).NE.NCJ) SQCDF = SQCDF+HOLD(I,J)*HOLD(I,J)
  100    CONTINUE
  140 CONTINUE
C
      SQCDF = SQRT(SQCDF)
      IF(OUT) WRITE (IW,9008)
      IF(OUT) CALL PREV(HOLD,E,NORB,L0,L1)
      RETURN
C
 9008 FORMAT(/5X,'CURRENT ORBITALS IN TERMS OF OLD ORBITALS')
      END
C*MODULE GVB     *DECK GVBDIA
      SUBROUTINE GVBDIA(A,VEC,EIG,SCR,IWRK,L0,L1,L2,LDV)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      DIMENSION A(L2),VEC(LDV,L1),EIG(L1),SCR(L1,8),IWRK(L1)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- DIAGONALIZE GVB FOCK MATRIX -A- WITHIN EACH SHELL -----
C
      CALL VCLR(VEC,1,LDV*L1)
      CALL VCLR(EIG,1,L1)
C
      IEND = 0
      DO 190 IHAM=1,NHAM+1
C
C            LOCATE STARTING/ENDING MO INDEX FOR THIS BLOCK
C
         ISTA = IEND+1
         IF(ISTA.GT.L0) GO TO 190
         IEND = 0
         IF(IHAM.LE.NHAM) THEN
            DO 110 I=NORB,1,-1
               IF(IEND.EQ.0  .AND.  NCONF(I).EQ.IHAM) IEND=I
  110       CONTINUE
         ELSE
            IEND = L0
         END IF
         NBLOCK = IEND + 1 - ISTA
C
C            COPY BLOCK FOR SHELL -IHAM- TO BEGINNING OF -A-
C
         DO 140 I=1,NBLOCK
            DO 130 J=1,I
               A(IA(I)+J)=A(IA(I+ISTA-1)+(J+ISTA-1))
  130       CONTINUE
  140    CONTINUE
C
C            DIAGONALIZE IT
C
         CALL GLDIAG(LDV,NBLOCK,NBLOCK,A,SCR,EIG(ISTA),
     *               VEC(1,ISTA),IERR,IWRK)
         IF(IERR.NE.0) THEN
            IF(MASWRK) WRITE(IW,*) 'FAILURE IN -GVBDIA-, IHAM=',IHAM
            CALL ABRT
            STOP
         END IF
C
C            MOVE THE ROTATION BLOCK ONTO THE DIAGONAL BLOCK
C
         IF(IHAM.EQ.1) GO TO 190
C
         DO 160 J=1,NBLOCK
            DO 150 I=1,NBLOCK
               VEC(NBLOCK-I+ISTA,J+ISTA-1)=VEC(NBLOCK-I+1,J+ISTA-1)
               VEC(NBLOCK-I+1   ,J+ISTA-1)=ZERO
  150       CONTINUE
  160    CONTINUE
  190 CONTINUE
      RETURN
      END
C*MODULE GVB     *DECK GVBFCK
      SUBROUTINE GVBFCK(JKUNIT,H,H2,S,TRANS,SCR,ERRMAT,WRK1,WRK2,DMP,
     *                  DODIIS,L1,L2,L3,OUT,SOSCF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL DODIIS,OUT,SOSCF,MFRZ
C
      DIMENSION H(L2),H2(L2),S(L2),TRANS(L3),SCR(L1),ERRMAT(L3),
     *          WRK1(L3),WRK2(L3)
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /MFRPAR/ MFRZ,NUMFRZ,IFRZ(MXAO)
C
      PARAMETER (ZERO=0.0D+00)
C
C     THIS ROUTINE MAKES THE GVB-PP FOCK MATRIX IN THE MO BASIS.
C     MODIFICATION OF J.MONTGOMERY ROUTINE FOR SINGLE-PAIR GVB.
C
C     THE COMBINED FOCK MATRIX IS RETURNED IN -H-
C     H(JI)=<J|H(J)-H(I)|I>, WHERE H(J),H(I) - FOCK OPERATORS
C           FOR SHELLS J AND I, ORBITALS J AND I ARE IN SHELLS J AND I
C     H(JI)=<J|H(J)|I> IF I IS UNOCCUPIED ORBITAL, J IS IN SHELL J
C     H(JI)=<J|H(J)|I> IF ORBITALS I AND J ARE FROM SAME SHELL J
C     H(JI)=<J|H(CORE)|I> IF BOTH I AND J ARE UNOCCUPIED
C
C        IF DIIS IS ACTIVE (DIIS WORKS ONLY FOR NPAIR=0 OR 1 SO FAR):
C     -H- WILL BE THE SINGLE FOCK OPERATOR IN AO BASIS
C     -ERRMAT- THE ERROR MATRIX IN THE AO BASIS
C
C      WHEN DIIS IS ACTIVE, BLOCKS PROPORTIONAL TO A**2 OR B**2
C      HAVE THESE FACTORS DIVIDED OUT TO AVOID POSSIBLE
C      PROBLEMS OF SLOW CONVERGENCE WHEN THEY BECOME SMALL.
C
C      ON ENTRY, -S- AND -TRANS- ARE OVERLAPS AND LCAO EXPANSION COEFS.
C
      DO 170 I = NORB+1,L1
         NCONF(I) = NCONF(NORB)+1
  170 CONTINUE
C
C        SKIP OVER CORE FOCK AND OPEN SHELL J/K OPERATORS
C
      LEN = L2
C
C     ----- READ CLOSED SHELL FOCK OPERATOR IN MO BASIS -----
C
      IF(NCO.GT.0) THEN
         CALL SEQREW(JKUNIT)
         CALL SEQADV(JKUNIT)
         NHAM1 = NHAM
         IF(NCO.GT.0) NHAM1 = NHAM - 1
         DO 180 I = 1,NHAM1
            CALL SEQADV(JKUNIT)
            CALL SEQADV(JKUNIT)
  180    CONTINUE
         CALL SQREAD(JKUNIT,H,LEN)
         IF(LEN.EQ.0) GO TO 800
C
C     ---- CONSTRUCT A QUASI-TOTAL FOCK OPERATOR ----
C     THE CODE HERE PRODUCES F = H + SUM 2*F(I)*J(I) - F(I)*K(I)
C     THIS SETS THE VIRTUAL-VIRTUAL BLOCK TO SOMETHING REASONABLE,
C     AND IT IS WHAT -GVBHAM- WOULD BE DOING TO SET UP H(CORE).
C
      ELSE
         CALL SEQREW(JKUNIT)
         CALL SQREAD(JKUNIT,H,L2)
         DO 220 IHAM=1,NHAM
            FACTA = F(IHAM)+F(IHAM)
            CALL SQREAD(JKUNIT,H2,L2)
            CALL DAXPY(L2,FACTA,H2,1,H,1)
            FACTB = -F(IHAM)
            CALL SQREAD(JKUNIT,H2,L2)
            CALL DAXPY(L2,FACTB,H2,1,H,1)
  220    CONTINUE
      END IF
C
      IF(OUT) THEN
         WRITE(IW,*) 'CORE FOCK MATRIX'
         CALL PRTRI(H,L1)
      END IF
C
C     --- READ SEPARATE FOCK MATRICES AND FORM SINGLE FOCK MATRIX ---
C
      NH1=2
      IF(NCO.EQ.0) NH1 = 1
      DO 390 IHAM=NH1,NHAM
C
         CALL SQREAD(JKUNIT,H2,LEN)
         IF(LEN.EQ.0) GO TO 800
C
         CALL DSCAL(L2,F(IHAM),H2,1)
         IF(OUT) THEN
            WRITE(IW,*) 'FOCK MATRIX NO.',IHAM
            CALL PRTRI(H2,L1)
         END IF
C
         DO 380 I=NCO+1,NORB
            IF (NCONF(I).NE.IHAM) GO TO 380
            DO 370 J=1,L1
               JI=J*(J-1)/2+I
               IF (J.LT.I) JI=I*(I-1)/2+J
               IF (.NOT.SOSCF) THEN
                  F1=F(IHAM)
                  F2=F(NCONF(J))-F(NCONF(I))
                  IF (DMP .GT. ZERO) THEN
                     F1 = F1+DMP
                     F2 = F2+DMP
                  END IF
               END IF
               IF (NCONF(J).LT.IHAM) THEN
                  IF (.NOT.SOSCF .AND. J.LE.NCO) THEN
                     H(JI)=(H(JI)-H2(JI))/F2
                  ELSE
                     H(JI)=H(JI)-H2(JI)
                  END IF
               ELSE
                  IF (.NOT.SOSCF .AND. J.GT.NORB) THEN
                     H(JI)=H2(JI)/F1
                  ELSE
                     H(JI)=H2(JI)
                  END IF
               END IF
  370       CONTINUE
  380    CONTINUE
  390 CONTINUE
C
      IF(OUT) THEN
         WRITE(IW,*) 'COMBINED FOCK MATRIX IN MO BASIS'
         CALL PRTRI(H,L1)
      END IF
C
C     IF MFRZ, ZERO OUT THE NON-DIAGONAL FROZEN ORBITAL ELEMENTS
C
      IF(MFRZ) CALL FRFOCK(H,L1)
      IF(MFRZ.AND.OUT) THEN
         WRITE(IW,*)'FOCK MATRIX IN MO BASIS AFTER ZEROING',
     *              ' THE FROZEN THE OFF-DIAG ELEMENTS'
         CALL PRTRI(H,L1)
      END IF
C
      IF(.NOT.DODIIS) RETURN
C
C        IF DIIS IS ACTIVE, -H- WILL BE THE SINGLE FOCK OPERATOR
C            AND -ERRMAT- THE ERROR MATRIX, BOTH IN THE AO BASIS.
C
C     ----- BACK-TRANSFORM GVB FOCK MATRIX TO AO BASIS -----
C     FAO = (SC) * HMO * (SC)-DAGGER
C
      CALL MTARBR(S,L1,TRANS,L1,WRK1,L1,1)
      CALL TRPOSQ(WRK1,L1)
C
      CALL DCOPY(L2,H,1,H2,1)
      CALL TFTRI(H,H2,WRK1,SCR,L1,L1,L1)
      IF(OUT) THEN
         WRITE(IW,*) 'COMBINED FOCK MATRIX IN AO BASIS'
         CALL PRTRI(H,L1)
      END IF
C
C     ----- MAKE ERROR MATRIX -----
C     IN THE MO BASIS, THIS IS ALL INTER-SHELL ELEMENTS OF FMO
C
      IJ = 0
      DO 260 I = 1,L1
         ISH = NCONF(I)
         DO 250 J = 1,I
            IJ = IJ + 1
            JSH = NCONF(J)
            IF (JSH.EQ.ISH) H2(IJ) = ZERO
  250    CONTINUE
  260 CONTINUE
C
C     ----- BACK TRANSFORM ERROR MATRIX TO AO BASIS -----
C     EXPAND FMO INTO AN ANTISYMMETRIC MATRIX, THEN BACKTRANSFORM.
C          AO ERR = (SC) * ERR * (SC)-DAGGER
C     NOTE THAT FROM ABOVE -WRK1- CONTAINS (SC)-DAGGER
C
      CALL EXPND(H2,ERRMAT,L1,1)
      CALL MRARBR(ERRMAT,L1,L1,L1,WRK1,L1,L1,WRK2,L1)
      CALL TRPOSQ(WRK1,L1)
      CALL MRARBR(WRK1,L1,L1,L1,WRK2,L1,L1,ERRMAT,L1)
C
      IF(OUT) THEN
         WRITE(IW,*) 'GVB-DIIS ERROR MATRIX IN AO BASIS'
         CALL PRSQ(ERRMAT,L1,L1,L1)
      END IF
      RETURN
C
  800 CONTINUE
      WRITE(IW,*) 'UNEXPECTED END OF FILE ON UNIT -JKUNIT- ... BOMBING'
      CALL ABRT
      STOP
      END
C*MODULE GVB     *DECK GVBGRD
      SUBROUTINE GVBGRD(H,GRAD,L1,L2,NPR,ORBGRD)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      DIMENSION H(L2),GRAD(NPR)
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- FORM ORBITAL GRADIENT FOR GVB-SOSCF CALCULATION FROM
C           MATRIX ELEMENTS OF COMBINED FOCK MATRIX IN MO BASIS
C           GRAD(JI) = (J | H(J)-H(I) | I>
C
C          H      - SINGLE COMBINED FOCK MATRIX IN MO BASIS
C          L1     - NUMBER OF BASIS FUNCTIONS
C          NPR    - NUMBER OF ROTATION PARAMETERS
C          ORBGRD - MAXIMUM ORBITAL GRADIENT COMPONENT
C
      DO 210 I = NORB+1,L1
         NCONF(I) = NCONF(NORB)+1
  210 CONTINUE
C
      K=0
      ORBGRD=ZERO
      DO 290 IHAM = 1,NHAM
         DO 280 J = 1,NORB
            IF (NCONF(J) .NE. IHAM) GO TO 280
            DO 270 I = NCO+1,L1
               IF (NCONF(I) .LE. IHAM) GO TO 270
               IJ=I*(I-1)/2+J
               K=K+1
               GRAD(K)=-H(IJ)
               IF (ABS(GRAD(K)).LT.1.0D-08) GRAD(K)=ZERO
  270       CONTINUE
  280    CONTINUE
  290 CONTINUE
      IMAX   = IDAMAX(NPR,GRAD,1)
      ORBGRD = ABS(GRAD(IMAX))
      RETURN
      END
C*MODULE GVB     *DECK GVBHAM
      SUBROUTINE GVBHAM(JKUNIT,DJ,DK,H,IA,L1,L2,OUT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL OUT
C
      DIMENSION DJ(L2),DK(L2),H(L2),IA(L2)
C
      COMMON /CONV  / DENTOL,EN,ETOT,EHF,EHFO,DIFF,ITER,ICALP,ICBET
      COMMON /GVBENG/ EK,VIR,TKIN(25),TPLUSV(25),ASYM
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
C     ----- CALCULATE A FOCK MATRIX FOR EACH FOCK OPERATOR -----
C     ----- CALCULATE T AND T + V
C
      EK = ZERO
      EHF = ZERO
      DO 100 K = 1,NHAM
         EK  = EK +TKIN(K)  *F(K)
         EHF = EHF+TPLUSV(K)*F(K)
  100 CONTINUE
C
C     ----- LOOP OVER FOCK OPERATORS -----
C
      NCORE = 1
      IF(NCO.EQ.0) NCORE=0
      DO 240 IHAM = 1,NHAM
         FACT = ONE/F(IHAM)
         CALL SEQREW(JKUNIT)
         CALL SQREAD(JKUNIT,H,L2)
         IF (OUT) THEN
            WRITE(IW,*) 'INCOMPLETE CORE FOCK OP. IN MO BASIS'
            CALL PRTRI(H,L1)
         END IF
         NHAM1 = NHAM
         IF (NCO.NE.0) NHAM1 = NHAM - 1
         IF (NHAM1.LT.1) GO TO 160
         DO 140 K = 1,NHAM1
            NX = K + NCORE
            NXIHAM = IA(NX)+IHAM
            IF (IHAM.GT.NX) NXIHAM = IA(IHAM)+NX
            CALL SQREAD(JKUNIT,DJ,L2)
            CALL SQREAD(JKUNIT,DK,L2)
            FACTA = ALPHA(NXIHAM)*FACT
            FACTB =  BETA(NXIHAM)*FACT
            IF (OUT) THEN
               WRITE(IW,9801) IHAM,K,FACT,FACTA,FACTB
 9801          FORMAT(' IHAM =',I3,' K =',I3,' FACT =',G12.5,
     *                ' FACTA =',G13.5,' FACTB =',G13.5)
               WRITE(IW,*) 'OPEN SHELL J OPERATOR FOR SHELL',K
               CALL PRTRI(DJ,L1)
               WRITE(IW,*) 'OPEN SHELL K OPERATOR FOR SHELL',K
               CALL PRTRI(DK,L1)
            END IF
            CALL DAXPY(L2,FACTA,DJ,1,H,1)
            CALL DAXPY(L2,FACTB,DK,1,H,1)
  140    CONTINUE
         IF (OUT) THEN
            WRITE(IW,*) 'TOTAL FOCK MATRIX FOR SHELL IHAM',IHAM
            CALL PRTRI(H,L1)
         END IF
C
C        ----- ADD THE TWO ELECTRON CONTRIBUTION TO TOTAL ENERGY -----
C
  160    CONTINUE
         DO 180 I = 1,NORB
            IF (NCONF(I) .NE. IHAM) GO TO 180
            IJ = IA(I)+I
            EHF = EHF+F(IHAM)*H(IJ)
  180    CONTINUE
C
C     ----- APPEND THE FOCK OPERATOR TO THE JK MATRIX FILE ----
C
         IDO = IHAM-1
         IF(IDO.GT.0) THEN
            DO 200 I = 1,IDO
               CALL SEQADV(JKUNIT)
  200       CONTINUE
         END IF
         CALL SQWRIT(JKUNIT,H,L2)
  240 CONTINUE
      RETURN
      END
C*MODULE GVB     *DECK GVBHES
      SUBROUTINE GVBHES(JKUNIT,HESS,HI,DIAGH,DIAGJ,DIAGK,
     *                  L1,L2,NPR,NHAMX)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      DIMENSION HI(L2),HESS(NPR),DIAGH(NHAMX,L1),DIAGJ(NHAMX,L1),
     *          DIAGK(NHAMX,L1)
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
      PARAMETER (TWO=2.0D+00, ONE=1.0D+00)
C
C     ----- CALCULATE DIAGONAL ELEMENTS OF INITIAL ORBITAL HESSIAN
C           FOR GVB-SOSCF
C
C         L1  = NUMBER OF BASIS FUNCTIONS
C         NPR = NUMBER OF ROTATION PARAMETERS
C         DIAGONAL ELEMENTS OF FOCK MATRICES ARE STORED IN DIAGH,
C         DIAGONAL ELEMENTS OF J AND K MATRICES IN DIAGJ AND DIAGK
C
      DO 100 I = NORB+1,L1
         NCONF(I) = NCONF(NORB)+1
  100 CONTINUE
C
      CALL VCLR(DIAGH,1,L1*NHAM)
      CALL VCLR(DIAGJ,1,L1*NHAM)
      CALL VCLR(DIAGK,1,L1*NHAM)
C
C     READ DIAGONAL ELEMENTS OF J AND K MATRICES
C
      NHAM1 = NHAM
      IF(NCO.GT.0) NHAM1 = NHAM - 1
      CALL SEQREW(JKUNIT)
      CALL SEQADV(JKUNIT)
      DO 290 IHAM = 1,NHAM1
         CALL SQREAD(JKUNIT,HI,L2)
         DO 230 IORB=1,L1
            IORBI=IORB*(IORB+1)/2
            DIAGJ(IHAM,IORB)=HI(IORBI)
  230    CONTINUE
         CALL SQREAD(JKUNIT,HI,L2)
         DO 240 IORB=1,L1
            IORBI=IORB*(IORB+1)/2
            DIAGK(IHAM,IORB)=HI(IORBI)
  240    CONTINUE
  290 CONTINUE
C
C     READ DIAGONAL ELEMENTS OF FOCK MATRICES
C
      DO 390 IHAM = 1,NHAM
         CALL SQREAD(JKUNIT,HI,L2)
         CALL DSCAL(L2,F(IHAM),HI,1)
         DO 370 IORB=1,L1
            IORBI=IORB*(IORB+1)/2
            DIAGH(IHAM,IORB)=HI(IORBI)
  370    CONTINUE
  390 CONTINUE
C
C     CALCULATE HESSIAN ELEMENTS
C     THE DIAGONAL HESSIAN ELEMENTS ARE EXACT EXCEPT FOR
C     CORE-VIRT ROTATIONS, SINCE TWO-ELECTRON INTEGRALS J AND K
C     ARE AVAILABLE ONLY FOR OPEN SHELLS
C
      K=0
      DO 490 IHAM = 1,NHAM
         DO 480 J = 1,NORB
            IF (NCONF(J) .NE. IHAM) GO TO 480
            DO 470 I = NCO+1,L1
               IF (NCONF(I) .LE. IHAM) GO TO 470
               K=K+1
               KII=NCONF(I)*(NCONF(I)+1)/2
               KJJ=NCONF(J)*(NCONF(J)+1)/2
               KIJ=NCONF(I)*(NCONF(I)-1)/2 + NCONF(J)
               IF (I.LE.NORB) THEN
                  HESS(K)=DIAGH(NCONF(J),I)-DIAGH(NCONF(I),I)
     *                -DIAGH(NCONF(J),J)+DIAGH(NCONF(I),J)
                  IF (NPAIR.GT.0) THEN
                  GAM=TWO*(ALPHA(KII)+ALPHA(KJJ)-TWO*ALPHA(KIJ))*
     *                DIAGK(NCONF(I)-1,J) + (BETA(KII)+BETA(KJJ)-
     *                TWO*BETA(KIJ))*(DIAGJ(NCONF(I)-1,J)+
     *                DIAGK(NCONF(I)-1,J))
                  HESS(K)=HESS(K)+GAM
                  END IF
               ELSE
                  HESS(K)= DIAGH(NCONF(J),I)-DIAGH(NCONF(J),J)
                  IF (NCONF(J).GT.1 .AND. NPAIR.GT.0) THEN
                     GAM=TWO*ALPHA(KJJ)*DIAGK(NCONF(J)-1,I)
     *                      + BETA(KJJ)*(DIAGJ(NCONF(J)-1,I)+
     *                                    DIAGK(NCONF(J)-1,I))
                     HESS(K)=HESS(K)+GAM
                  END IF
               END IF
  470       CONTINUE
  480    CONTINUE
  490 CONTINUE
C
C     INVERT THE DIAGONAL HESSIAN
C
      HSSTOL = 0.005D+00
      DO 500 IPR=1,NPR
         IF(ABS(HESS(IPR)).GT.HSSTOL) THEN
            HESS(IPR) = ONE/HESS(IPR)
         ELSE
            HESS(IPR) = ONE/HSSTOL
         END IF
  500 CONTINUE
C
      RETURN
      END
C*MODULE GVB     *DECK GVBITR
      SUBROUTINE GVBITR(TRANS,DJ,DK,H,H2,E,EIG,XP,XK,IXPK,NIJS,NKLS,
     *                  S,Q,SCR,IWRK,IA,DSH,XINTS,GHONDO,DDIJ,
     *                  MXFOCK,NINTMX,L0,L1,L2,L3,NSH2,MAXG,INTTYP,
     *                  DIRSCF,SCHWRZ,WRK1,WRK2,WRK3,ADIIS,XDIIS,
     *                  PDIIS,BDIIS,IODIIS,HESS,GRAD,PGRAD,DISPL,
     *                  WRK4,WRK5,WRK6,DIAGH,DIAGJ,DIAGK,NHAMX,
     *                  MAXIT1,MAXIT2,MXDII3,MAXIO,NPR)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXAO=8192,
     *           MXCAV=5810, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      LOGICAL OUT,SOME,NOTOPN,GOPARR,DSKWRK,MASWRK,SVDSKW
      LOGICAL CVGED,CVGING,CVDENS,CVENGY,CVDIIS
      LOGICAL CVGSVP,CVGSVPDN,CVGSVPDN1,CVGSVPNR
      LOGICAL EXTRAH,DAMPH,VSHIFT,RSTRCT,DODIIS,ROTDMP
      LOGICAL DIRSCF,FDIFF,SCHWRZ,MINMEM,SOSCF,SOSTART,FT15OP
      LOGICAL IPCFP, IZRFLS
C
      DIMENSION XP(NINTMX),XK(NINTMX),IXPK(NINTMX),NIJS(NINTMX),
     *          NKLS(NINTMX),DSH(NSH2),XINTS(NSH2),GHONDO(MAXG),
     *          EIG(L1),IA(L1),SCR(L1,8),E(L1),TRANS(L1,L1),IWRK(L1),
     *          WRK1(L3),WRK2(L3),WRK3(L3),S(L2),Q(L3),DJ(*),DK(*),
     *          H(L2),H2(L2),ADIIS(MXDII3),XDIIS(MAXIT1),PDIIS(MAXIT1),
     *          BDIIS(MAXIT2),IODIIS(MAXIO),DDIJ(*),
     *          HESS(NPR),GRAD(NPR),PGRAD(NPR),DISPL(NPR),
     *          WRK4(NPR),WRK5(NPR),WRK6(NPR),
     *          DIAGH(NHAMX,L1),DIAGJ(NHAMX,L1),DIAGK(NHAMX,L1)
C
      COMMON /ACONV / RRSHFT,EXTTOL,DMPTOL,VSHTOL,IEXTIN
      COMMON /ATHEAT/ ATHEAT
      COMMON /CONV  / DENTOL,EN,ETOT,EHF,EHF0,DIFF,ITER,ICALP,ICBET
      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
      COMMON /DMPING/ SHIFTO,SHIFTV,DMPCUT,SWDIIS,DIRTHR
      COMMON /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      COMMON /GVBCI / CIHAM(91),CILOW(12),JPAIR,KONE,KTWO,KCORB(2,12)
      COMMON /GVBENG/ EK,VIR,TKIN(25),TPLUSV(25),ASYM
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /FMCOM / X(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FUNCT / ENERGY,EGRAD(3,MXATM)
      COMMON /INDDIP/ TMUX,TMUY,TMUZ,TMUXD,TMUYD,TMUZD,MINMEM
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /MOLMEC/ HTYPE(4),NHCO(4,200),NNHCO,ITYPE
      COMMON /MPCGEO/ GEO(3,MXATM)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FEX,UNZOR,PORT
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMEFP/ PB_MUL,PC_MUL,PB_POL,PC_POL,P_FF,P_NF,P_FN,P_I,
     *                P_IBIS,P_J,P_NUCC,P_NUCCBIS,ENPCM
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IRPCM,IEF,IP_F,NFMOPCM
      COMMON /PCMPNT/ LXYZRE,LSSFE,LLIST,LALPCM,LRINPCM,LINAPCM,LINFPCM,
     *                LINIPCM,LMEPCM,LAXYZCT,LXYZCT2,LNVERT,LQSN,LQSE,
     *                LQFS,LQIND,LISPHE,LVAD,LQOR,LVECMUL,LAIPRJ,LFIPRJ
     *                ,LPEL
      COMMON /PCMPRT/ GCAVP,GCAVS,GDISP,GREP,EHFGAS
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /SCFWFN/ AROHF(3),BROHF(3),PACAVO(6),IACAVO,IUHFNO,MVOQ
      COMMON /SVPSCM/ IADDR_SVP,ILAST_SVP
      COMMON /SVPNRG/ ATMSVP,ELCSVP,SVPSVP,SURSVP,VOLPLM,VOLSVP,RXNSVP
      COMMON /SVPOPT/ ISVP,NSVP
      COMMON /SVPTGR/ NVLPL,MTHSVP,IQORP,NSTAGE,IOPPRD,IRENRM,ISHAPE,
     *                ITRCAV,NDRCAV,ICVICK,IFRQUP,IRDRF,IONSTR,
     *                INTCAV,NPTLEB,NPTTHE,NPTPHI,IROTGR,ITRNGR,
     *                ICRD,LINEQ,IPNRF,IPRRF,IPNDN,IPRDN,IEX
      COMMON /SVPPAR/ CAVC(3,MXCAV),CAVW(MXCAV),CAVQ(MXCAV),
     1                CAVN(3,MXCAV),SPHW(MXCAV),CAVP(MXCAV),
     2                CENCAV(3),RNTCAV(3),NCAV,IQP
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, TEN=1.0D+01)
      PARAMETER (HALF=0.5D+00, PT2=0.2D+00, TWOPT2=2.2D+00)
      PARAMETER (TOLZ=1.0D-08, TOLE=1.0D-05, FIVE=5.0D+00)
C
      SAVE FT15OP
C
      DATA FT15OP/.FALSE./
      DATA DEBUG,DBGME,CHECK/8HDEBUG   ,8HGVB     ,8HCHECK   /
      DATA GVB/8HGVB     /
      DATA OPTMZE/8HOPTIMIZE/
      DATA DMPTLC/1.0D-02/
      DATA HESSN/8HHESSIAN /
      DATA NONE/4HNONE/
      DATA ENER/8HENERGY  /
C
      IZRFLS = IZRF.EQ.1
      IF(IPCM.EQ.1) EN = EN - UNZ
C
C     ----- ROHF-GVB ITERATION DRIVER -----
C
      OUT = NPRINT.EQ.5  .OR.  EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBGME
     *      .AND. MASWRK
      SOME = NPRINT.NE.-5 .AND. MASWRK
      NPASS = 0
C
C     ----- ALLOCATE A FILE FOR SAVING THE J AND K MATRICES -----
C
      SVDSKW = DSKWRK
      IF (.NOT.DIRSCF) DSKWRK = .TRUE.
      JKUNIT = 23
      CALL SEQOPN(JKUNIT,'JKFILE','UNKNOWN',.FALSE.,'UNFORMATTED')
      DSKWRK = SVDSKW
C
C     ----- START THE CLOCK -----
C
      TIME0 = ZERO
      TIME1 = ZERO
      CALL TSECND(TIME0)
      TLEFTS = TIMLIM - TIME0
C
C     ----- INITIALIZE SECTION TIMERS -----
C
      TIMITT= ZERO
      TIMEF = ZERO
      TIMEG = ZERO
      TIMEM = ZERO
      TIMEO = ZERO
      TIMDUM= ZERO
C
C     ---- INITIALIZE SOME CONSTANTS FOR CONVERGENCE CONTROL -----
C     NOTE THAT DIIS IS ACTIVE FOR NPAIR=0 AND 1 ONLY.
C
      CVGING = .FALSE.
      CVGED  = .FALSE.
C
C     CVGSVP KEEPS TRACK OF WHETHER SVP CHARGES/DIPOLES APPEAR CONVERGED
C
      CVGSVP    = .FALSE.
      CVGSVPDN  = .FALSE.
      CVGSVPDN1 = .FALSE.
      CVGSVPNR  = .FALSE.
C
      IF(DIRSCF) THEN
         DENTOL = TWO * CONVHF
         ENGTOL = 1.0D-08
         DIITOL = 1.0D-06
      ELSE
         DENTOL = CONVHF
         ENGTOL = 1.0D-09
         DIITOL = 1.0D-07
      END IF
      IF(RUNTYP.EQ.OPTMZE  .AND.  NEVALS.GT.0) THEN
         GRMS = DDOT(3*NAT,EGRAD,1,EGRAD,1)
         GRMS = SQRT(GRMS/(3*NAT))
                               CFACT =  1.0D+00
         IF(GRMS.GT.0.005D+00) CFACT =  5.0D+00
         IF(GRMS.GT.0.020D+00) CFACT = 20.0D+00
         IF(GRMS.GT.0.100D+00) CFACT = 50.0D+00
         DENTOL = CFACT*DENTOL
         ENGTOL = CFACT*ENGTOL
         DIITOL = CFACT*ENGTOL
      END IF
      DTOL2 = TWO * DENTOL
      ETOL2 = TEN * ENGTOL
      CVDIIS = .FALSE.
C
      EXTRAH = MOD(MCONV, 2) .EQ.  1
      DAMPH  = MOD(MCONV, 4) .GE.  2
      VSHIFT = MOD(MCONV, 8) .GE.  4
      RSTRCT = MOD(MCONV,16) .GE.  8
      DODIIS = MOD(MCONV,32) .GE. 16
      SOSCF  = MOD(MCONV,128).GE. 64
      IF(NPAIR.GT.1) DODIIS=.FALSE.
      ROTDMP = DAMPH
C
      DO 10 I=1,NSETO
         IF(NO(I).EQ.1) GO TO 10
            RSTRCT = .TRUE.
            IF(SOME) WRITE(IW,9010)
   10 CONTINUE
C
      DAMP = ZERO
      DAMP0 = ZERO
      IF (DAMPH .OR. VSHIFT) DAMP = MAX(ONE,DMPCUT)
      IF (DMPCUT .LE. ZERO) DMPCUT = ZERO
      ITERV = 0
      ITER = 0
      ICALP = 0
      RRSHFT = ZERO
      DIFF = ZERO
      DIFF0 = DIFF
      DIFFP = ZERO
      DELE  = ZERO
      DELE0 = ZERO
      DEAVG = ZERO
      EHF  = ZERO
      EHF0 = ZERO
      ETOT = ZERO
C
C        MOPAC HAS A MOLECULAR MECHANICS PEPTIDE BOND CORRECTION
C
      EHNCO = ZERO
      IF(MPCTYP.NE.NONE) THEN
         DO I=1,NNHCO
            CALL DIHED(GEO,NHCO(1,I),NHCO(2,I),
     *                     NHCO(3,I),NHCO(4,I),ANGLE,SIGNA)
            SINANG = SIN(ANGLE)
            EHNCO = EHNCO + HTYPE(ITYPE)*SINANG*SINANG
         ENDDO
         EHNCO = EHNCO/627.52791D+00
      END IF
C
      CALL VCLR(CIHAM,1,91)
C
      ERDIIS = ZERO
      ITDIIS = 1
      DMPDII = 0.02D+00
      NOTOPN = .TRUE.
C
C        RUNS W/O PAIRS SEEM TO BENEFIT FROM ONE CONVENTIONAL
C        MATRIX DIAGONALIZATION STEP.  SO TOO WILL HESSIAN
C        RUNS, AT ALL NUMERIC DISPLACEMENTS, SINCE THE DIAGS
C        WILL BREAK THE ORBITAL SYMMETRY FROM VIB 0 UP.  IT
C        SEEMS BEST TO START SOSCF AT ONCE FOR ALL PAIR RUNS.
C
      ITSO=0
      NSO=1
      IF(NEVALS.GT.0) NSO=0
      IF(RUNTYP.EQ.HESSN) NSO=1
      IF(NPAIR.EQ.1) NSO=0
      ORBGRD=ZERO
C
C     ----- INITIALIZE DISK FOR EXTRAPOLATION -----
C
      NDAF = 71
      LNOC = L1*NORB
      IF(LNOC.LT.L2) LNOC = L2
      CALL WDISK(TRANS,TRANS,TRANS,NDAF,LNOC)
C
C     ----- PRINT ITERATION HEADER -----
C
      IF(SOME) THEN
         WRITE(IW,9018) MAXIT,NPUNCH,DENTOL,EN,
     *                  EXTRAH,DAMPH,VSHIFT,RSTRCT,DODIIS,SOSCF
         IF(SOSCF) WRITE(IW,9017) NPR,SOGTOL
      END IF
C
      IF(EXETYP.EQ.CHECK) GO TO 700
C
C     ----- READ OVERLAP INTEGRALS -----
C     ----- READ STARTING ORBITALS -----
C     ----- READ SYMMETRY ADAPTED -Q- MATRIX -----
C
      CALL DAREAD(IDAF,IODA,S    ,L2,12,0)
      CALL DAREAD(IDAF,IODA,TRANS,L3,15,0)
      CALL DAREAD(IDAF,IODA,Q    ,L3,45,0)
      CALL VCLR(E,1,L1)
C
C     ----- EXCHANGE INTEGRALS FOR DIRECT SCF THRESHOLD TESTS -----
C
      NINT=0
      NSCHWZ=0
      IF(DIRSCF) THEN
         FDIFF=.FALSE.
         IF(MASWRK) WRITE(IW,9030) SCHWRZ,FDIFF
         IF(SCHWRZ) CALL EXCHNG(XINTS,GHONDO,DDIJ,NSH2,MAXG,INTTYP)
      END IF
C
      NFT15=15
      IF(SOSCF) THEN
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         IF(.NOT.FT15OP) THEN
            FT15OP = .TRUE.
            CALL SEQOPN(NFT15,'WORK15','UNKNOWN',.FALSE.,'UNFORMATTED')
         ELSE
            CALL SEQREW(NFT15)
         END IF
         DSKWRK = SVDSKW
      END IF
C
C     -- SET UP COMBINED EFP/PCM
C
      IF(IPCM.EQ.1 .AND. IEFP.EQ.1) THEN
         IPCM=0
         IEFP=0
         IPCFP=.TRUE.
      ELSE
         IPCFP=.FALSE.
      END IF
C
C     -- MEMORY FOR ANY OF THE SOLVATION METHODS
C
      CALL VALFM(LOADFM)
      LAST   = LOADFM + 1
C
      IF (IZRFLS .OR. IEFP.EQ.1 .OR. IPCFP .OR.
     *    IPCM.EQ.1 .OR. NBUFMO.GT.0  .OR.  ISVP.EQ.1) THEN
                         LENW2= 0
         IF(NBUFMO.GT.0) LENW2=L2
         LWRK1  = LAST
         LWRK2  = LWRK1  + L2
         LAST   = LWRK2  + LENW2
      ELSE
         LWRK1  = LAST
         LWRK2  = LAST
      END IF
C
C        FOR PCM AND COMBINED PCM/EFP
C            THE ORDER MUST BE PCM, EFP, PCM/EFP
C
      CALL PCMMEM3(IPCFP,L2,LAST,NTSPAR,LSOL,LCHG,LDMATM,LQPOT,LQFLD,
     *             LVPOT,LSE,LDE,LQ0,LQ1,LQ2,LQ3,LD0,LQA,LDIMAT,LQREP,
     *             LTMP0,LTMP1,LTMP2,LIPVT,LPOTTMP,LRMUL,LCHG2,
     *             LASCCRD,LASCCHG,LASCDIP,LASCQAD,LDISV,LDIS1,LDIS2,
     *             LCQEF,LQEFF,LELD,LXDINT,LYDINT,LZDINT,LEFLD,MADD,
     *             LABFLD,LVECMP,LVEC_2,LVEC_S,LAPROJ,LFPROJ)
      LSCR   = LAST
      IF(IPCM.EQ.1.OR.IPCFP) THEN
         LCOL1  = LSCR   + L2
         LPOT   = LCOL1  + L1
         LBK    = LPOT   + L2
         LAST   = LBK    + L3
      ELSE
         LCOL1  = LAST
         LPOT   = LAST
         LBK    = LAST
      END IF
C
C         MEMORY FOR SVP SOLVENT MODEL
C
      IF(ISVP.EQ.1) THEN
         NOPER=2
         IF(IQP.EQ.0.AND.MTHSVP.EQ.1) NOPER=1
         LENSVP=10*NCAV+NOPER*NCAV**2
         LSVP=LAST
         LAST=LSVP+LENSVP
      ELSE
         LSVP=LAST
         LENSVP=0
      END IF
C          AND INITIALIZATIONS FOR -SVP- MODEL
      IUPRXN = -1
      IF(ISVP.EQ.1) THEN
         IADDR_SVP=LSVP
         ILAST_SVP=LSVP+LENSVP-1
         CALL SVPSET
      END IF
C
      IF (IZRFLS .OR. IEFP.EQ.1 .OR. IPCFP .OR.
     *    IPCM.EQ.1 .OR. NBUFMO.GT.0  .OR.  ISVP.EQ.1) THEN
         NEED = LAST-LOADFM-1
         IF(MASWRK) WRITE(IW,9020) NEED
         CALL GETFM(NEED)
      END IF
C
C     -- READ INVERSE D MATRIX FOR PCM
C
      IF((IPCM.EQ.1.OR.IPCFP).AND.IPCMIT.EQ.0) THEN
         NTS3 = NTS*NTS
         CALL SEQREW(NFT26)
         CALL DERCHK(NDER)
         IF(NDER.EQ.1.AND.IPCDER.EQ.1) THEN
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
         END IF
         IF(IEF.GT.0.AND.IEF.LT.3) THEN
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
            CALL SQREAD(NFT26,X(LSE),NTS3)
            CALL SQREAD(NFT26,X(LDE),NTS3)
         END IF
         CALL SQREAD(NFT26,X(LDMATM),NTS3)
         CALL SEQREW(NFT26)
      END IF
C
C     ----- SAVE PRISTINE ONE ELEC. HAM. FOR SOLVENT MODELS -----
C     ----- IF SCRF READ IN DIPOLE INTEGRALS -----
C
      IF(IZRFLS .OR. IEFP.EQ.1. OR. IPCM.EQ.1 .OR.
     *   IPCFP .OR. NBUFMO.GT.0  .OR.  ISVP.EQ.1) THEN
         CALL DAREAD(IDAF,IODA,H2,L2,11,0)
         IF(.NOT.(IPCFP.AND.RUNTYP.EQ.ENER)) THEN
            CALL DAWRIT(IDAF,IODA,H2,L2,87,0)
         END IF
      END IF
      IF(IZRFLS) THEN
         CALL DAREAD(IDAF,IODA,X(LXDINT),L2,95,0)
         CALL DAREAD(IDAF,IODA,X(LYDINT),L2,96,0)
         CALL DAREAD(IDAF,IODA,X(LZDINT),L2,97,0)
      END IF
C
C     ----- PRINT ITERATION HEADER -----
C
      IF (MASWRK) THEN
      IF(DAMPH  .OR.  VSHIFT) THEN
         IF(SOSCF) THEN
            IF(.NOT.DIRSCF) WRITE(IW,9040)
            IF(     DIRSCF) WRITE(IW,9050)
         ELSE
            IF(.NOT.DIRSCF) WRITE(IW,9045)
            IF(     DIRSCF) WRITE(IW,9055)
         END IF
      ELSE
         IF(SOSCF) THEN
            IF(.NOT.DIRSCF) WRITE(IW,9060)
            IF(     DIRSCF) WRITE(IW,9070)
         ELSE
            IF(.NOT.DIRSCF) WRITE(IW,9065)
            IF(     DIRSCF) WRITE(IW,9075)
         END IF
      END IF
      END IF
      CALL FLSHBF(IW)
C
C           ********************
C     ***** BEGIN GVB ITERATIONS *****
C           ********************
C
      CALL TSECND(TIMITT)
  100 CONTINUE
      IF(OUT) WRITE(IW,FMT='(1X,''***** ITER='',I4)') ITER
C
C
C     * * * * * SOLVATION MODEL CORRECTIONS TO 1E- TERM * * * * *
C
C                ----- ONSAGER REACTION FIELD -----
C     GET PRISTINE ONE ELEC. HAM., CALC. INDUCED DIPOLES, ADD PERTURB.
C     -H- HOLDS TOTAL DENSITY,
C     -H2- HOLDS CORE HAMILTONIAN PLUS SCRF PERTURBATION,
C     AND WILL BE SAVED ON THE DAF INSIDE RCNFLD IN RECORD 11.
C
      IF (IZRF.EQ.1) THEN
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         CALL DENDD1(H,H2,L2)
         CALL DAREAD(IDAF,IODA,H2,L2,87,0)
         CALL RCNFLD(H2,DUMMY,DUMMY,H,DUMMY,DUMMY,
     *               X(LXDINT),X(LYDINT),X(LZDINT),L2,L1,0)
         DSKWRK = SVDSKW
      END IF
C
C     ----- SURFACE AND VOLUME POLARIZATION FOR ELECTROSTATICS -----
C        -H- USED FOR TOTAL DENSITY, AND -H2- FOR H + SVP CORRECTIONS
C
      IF(ISVP.EQ.1) THEN
         CALL DENDD1(H,H2,L2)
         ITERSVP = ITER+1
         CALL DAREAD(IDAF,IODA,H2,L2,87,0)
         CALL SVPCHG(H,DUMMY,H2,X(LWRK1),IADDR_SVP,ILAST_SVP,
     1               L2,L1,ESVP,ITERSVP,
     2               CVGSVPDN,CVGSVPDN1,CVGSVPNR,IUPRXN,CHGDIF,CVGSVP)
         CALL DAWRIT(IDAF,IODA,H2,L2,11,0)
      END IF
C
C           ----- EFFECTIVE FRAGMENT POTENTIAL MODEL -----
C
C  -- IF FRAG GET CURRENT DENSITY MATRIX AND PRISTINE ONE --
C  -- ELEC. HAM., CALC. INDUCED DIPOLES, ADD PERTURBATION --
C     -H- HOLDS TOTAL DENSITY,
C     -H2- HOLDS CORE HAMILTONIAN PLUS IND. DIPOLE PERTURBATION,
C     AND WILL BE SAVED ON THE DAF INSIDE EFPCMP IN RECORD 11.
C
      IF(NBUFMO.GT.0) THEN
         CALL DAREAD(IDAF,IODA,H2,L2,87,0)
      END IF
C
      IF (IEFP.EQ.1) THEN
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         CALL DENDD1(H,H2,L2)
         CALL DAREAD(IDAF,IODA,H2,L2,87,0)
C
         IF(NBUFMO.GT.0) THEN
           CALL BUFDEN(H2,X(LWRK1),L1,L2)
           CALL VSUB(X(LWRK1),1,H,1,H,1,L2)
         END IF
C
         CALL EFPCMP(H2,DUMMY,DUMMY,H,DUMMY,DUMMY,
     *               X(LXDINT),X(LYDINT),X(LZDINT),X(LEFLD),
     *               X(MADD),X(LABFLD),L2,L1,0)
         DSKWRK = SVDSKW
      END IF
C
C     ----- POLARIZABLE CONTINUUM MODEL CALCULATION -----
C
C  -- IF PCM GET CURRENT DENSITY MATRIX AND PRISTINE ONE  --
C  -- ELEC. HAM., CALC. SURFACE CHARGES, ADD PERTURBATION --
C     -H- HOLDS TOTAL DENSITY,
C     -H2- HOLDS CORE HAMILTONIAN PLUS APPARENT SURFACE CHARGE
C     PERTURBATION, AND IS WRITTEN TO DAF 11 INSIDE PCMFLD
C
      IF(IPCM.EQ.1) THEN
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
C
C        -- ADD PCM 1-E INTERACTION MATRIX
C
         IF(ITER.EQ.0.AND.IPCMIT.EQ.0) THEN
            IF(IP_F.EQ.1) THEN
             CALL VALFM(LOADFM)
             LVEC    = LOADFM   + 1
             LVECI   = LVEC     + L2
             LVECTEST= LVECI    + L2
             LAST    = LVECTEST + L2
             NEEDEF = LAST - LOADFM -1
             CALL GETFM(NEEDEF)
             CALL ENUADD(X(LDMATM),X(LVEC),X(LVECI),X(LVECTEST),
     *                   L2,X(LAXYZCT),
     *                   X(LAXYZCT+MXTS),X(LAXYZCT+MXTS*2),
     *                   X(LAXYZCT+MXTS*3),X(LQSN),X(LQFS),X(LVECMUL))
             CALL DAWRIT(IDAF,IODA,X(LVECI),L2,263,0)
             CALL RETFM(NEEDEF)
            END IF
            EN = EN + UNZ
C
            CALL JMAT(X(LSOL),X(LSCR),X(LPOT),X(LDIS1),X(LDIS2),L2,
     *                X(LAXYZCT),X(LAXYZCT+MXTS),X(LAXYZCT+MXTS*2),
     *                X(LQSN),X(LQOR))
            CALL DAWRIT(IDAF,IODA,X(LSCR),L2,264,0)
         END IF
C
         CALL DENDD1(H,H2,L2)
         CALL DAREAD(IDAF,IODA,H2,L2,87,0)
C
         CALL DCOPY(L2,H,1,X(LWRK1),1)
         IF(NBUFMO.GT.0) THEN
           CALL BUFDEN(H2,X(LWRK1),L1,L2)
           CALL VSUB(X(LWRK1),1,H,1,X(LWRK1),1,L2)
         END IF
C
         CALL PCMFLD(H2,DUMMY,H,X(LSOL),X(LCHG),X(LSCR),X(LDISV),
     *               X(LDIS1),X(LDIS2),X(LBK),X(LCOL1),X(LDMATM),
     *               X(LCQEF),X(LQEFF),X(LELD),X(LQPOT),X(LQFLD),
     *               X(LVPOT),X(LSE),X(LDE),
     *               X(LQ0),X(LQ1),X(LQ2),X(LQ3),
     *               X(LD0),X(LQA),X(LDIMAT),X(LQREP),
     *               X(LTMP0),X(LTMP1),X(LTMP2),X(LIPVT),
     *               X(LPOTTMP),X(LRMUL),TCH,QET,QETN,QESC,1,NFT27,
     *               .FALSE.,L1,L2,NTSPAR,X(LAXYZCT),X(LAXYZCT+MXTS),
     *               X(LAXYZCT+MXTS*2),X(LXYZCT2),X(LXYZCT2+MXTS),
     *               X(LXYZCT2+MXTS*2),X(LAXYZCT+MXTS*3),
     *               X(LQSN),X(LQSE),X(LQFS),X(LVECMUL),
     *               X(LCHG2),X(LQIND+MXTS),X(LPEL),X(LASCCRD),
     *               X(LASCCHG),X(LASCDIP),X(LASCQAD),X(LXYZRE),
     *               X(LXYZRE+MXSP),X(LXYZRE+MXSP*2),X(LXYZRE+MXSP*3),
     *               X(LISPHE),X(LLIST))
         DSKWRK = SVDSKW
      END IF
C
C               ------ COMBINED PCM/EFP METHOD ------
C
C     -- IF EFP/PCM GET CURRENT DENSITY MATRIX AND PRISTINE    --
C     -- ONE ELEC. HAM., ADD PERTURBATION                      --
C     --  H   HOLDS TOTAL DENSITY,                             --
C     --  H2  HOLDS CORE HAMILTONIAN PLUS CHARGE AND DIPOLE    --
C     -- PERTURBATION, AND IS WRITTEN TO 11 IN PCMFLD & CHGIND --
C
      IF(IPCFP) THEN
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
C
C        -- ADD PCM 1-E INTERACTION MATRIX
C
         IF(ITER.EQ.0.AND.IPCMIT.EQ.0) THEN
            IF(IP_F.EQ.1) THEN
             CALL VALFM(LOADFM)
             LVEC    = LOADFM   + 1
             LVECI   = LVEC     + L2
             LVECTEST= LVECI    + L2
             LAST    = LVECTEST + L2
             NEEDEF = LAST - LOADFM -1
             CALL GETFM(NEEDEF)
             CALL ENUADD(X(LDMATM),X(LVEC),X(LVECI),X(LVECTEST),
     *                   L2,X(LAXYZCT),
     *                   X(LAXYZCT+MXTS),X(LAXYZCT+MXTS*2),
     *                   X(LAXYZCT+MXTS*3),X(LQSN),X(LQFS),X(LVECMUL))
             CALL DAWRIT(IDAF,IODA,X(LVECI),L2,263,0)
             CALL RETFM(NEEDEF)
            END IF
            EN = EN + UNZ
C
            CALL JMAT(X(LSOL),X(LSCR),X(LPOT),X(LDIS1),X(LDIS2),L2,
     *                X(LAXYZCT),X(LAXYZCT+MXTS),X(LAXYZCT+MXTS*2),
     *                X(LQSN),X(LQOR))
            CALL DAWRIT(IDAF,IODA,X(LSCR),L2,264,0)
         END IF
C
         CALL DENDD1(H,H2,L2)
         CALL DAREAD(IDAF,IODA,H2,L2,87,0)
C
         CALL DCOPY(L2,H,1,X(LWRK1),1)
         IF(NBUFMO.GT.0) THEN
           CALL BUFDEN(H2,X(LWRK1),L1,L2)
           CALL VSUB(X(LWRK1),1,H,1,X(LWRK1),1,L2)
         END IF
C
         CALL PCMFLD(H2,DUMMY,H,X(LSOL),X(LCHG),X(LSCR),X(LDISV),
     *               X(LDIS1),X(LDIS2),X(LBK),X(LCOL1),X(LDMATM),
     *               X(LCQEF),X(LQEFF),X(LELD),X(LQPOT),X(LQFLD),
     *               X(LVPOT),X(LSE),X(LDE),
     *               X(LQ0),X(LQ1),X(LQ2),X(LQ3),
     *               X(LD0),X(LQA),X(LDIMAT),X(LQREP),
     *               X(LTMP0),X(LTMP1),X(LTMP2),X(LIPVT),
     *               X(LPOTTMP),X(LRMUL),TCH,QET,QETN,QESC,1,NFT27,
     *               IPCFP,L1,L2,NTSPAR,X(LAXYZCT),X(LAXYZCT+MXTS),
     *               X(LAXYZCT+MXTS*2),X(LXYZCT2),X(LXYZCT2+MXTS),
     *               X(LXYZCT2+MXTS*2),X(LAXYZCT+MXTS*3),
     *               X(LQSN),X(LQSE),X(LQFS),X(LVECMUL),
     *               X(LCHG2),X(LQIND+MXTS),X(LPEL),X(LASCCRD),
     *               X(LASCCHG),X(LASCDIP),X(LASCQAD),X(LXYZRE),
     *               X(LXYZRE+MXSP),X(LXYZRE+MXSP*2),X(LXYZRE+MXSP*3),
     *               X(LISPHE),X(LLIST))
C
         CALL CHGIND(H2,DUMMY,DUMMY,X(LWRK1),DUMMY,DUMMY,
     *               X(LDMATM),X(LXDINT),X(LYDINT),X(LZDINT),
     *               X(LEFLD),X(MADD),X(LABFLD),X(LVECMP),X(LVEC_2),
     *               X(LVEC_S),L2,L1,0,
     *               X(LAXYZCT),X(LAXYZCT+MXTS),X(LAXYZCT+MXTS*2),
     *               X(LAXYZCT+MXTS*3),X(LQSN),X(LQSE),X(LQFS),
     *               X(LQIND),X(LAIPRJ),X(LFIPRJ))
         DSKWRK = SVDSKW
      END IF
C
C   SUBTRACT INTERACTION OF THE EFP WITH THE BUFFER ZONE DENSITY.
C
C   HERE WE DEFINE EEFP AS THE INTERACTION OF THE EFP(MULTIPOLE
C   POINTS AND POLARIZABLE POINTS) WITH THE BUFFER ZONE DENSITY
C   I.E. (1)BUFFER ELECTRON -- MULTIPOLE POINTS
C        (2)BUFFER ELECTRON -- POLARIZABLE POINTS(INDUCED DIPOLE)
C   ANY PCM ASC ASSOCIATED INTERACTION IS NOT CONSIDERED AS A
C   PART OF EEFP.
C
      IF(NBUFMO.GT.0) THEN
        SVDSKW = DSKWRK
        DSKWRK = .FALSE.
        CALL DAREAD(IDAF,IODA,X(LWRK1),L2,319,0)
        CALL VSUB(X(LWRK1),1,H2,1,X(LWRK1),1,L2)
        IF(IPCFP.OR.IPCM.EQ.1)
     *     CALL VSUB(X(LSOL),1,X(LWRK1),1,X(LWRK1),1,L2)
C
        IF(IPCMIT.EQ.0) THEN
           IF(IPCM.EQ.1.OR.IPCFP) THEN
              CALL DAREAD(IDAF,IODA,X(LSCR),L2,263,0)
              CALL VSUB(X(LSCR),1,X(LWRK1),1,X(LWRK1),1,L2)
              CALL DAREAD(IDAF,IODA,X(LSCR),L2,264,0)
              CALL VSUB(X(LSCR),1,X(LWRK1),1,X(LWRK1),1,L2)
           END IF
           IF(IPCFP) THEN
              CALL VSUB(X(LVEC_2),1,X(LWRK1),1,X(LWRK1),1,L2)
           END IF
        END IF
C
        CALL BUFDEN(H2,X(LWRK2),L1,L2)
        EEFP= TRACEP(X(LWRK2),X(LWRK1),L1)
        DSKWRK = SVDSKW
      ELSE
        EEFP = ZERO
      END IF
C
C     ----- DAMP AND EXTRAPOLATE THE EIGENVECTORS IF NECESSARY -----
C
C     NOTE THAT THIS IS DONE TO THE VECTORS, NOT A FOCK MATRIX,
C     BECAUSE THERE IS NEVER A SINGLE COMBINED FOCK MATRIX.
C     DAMPD COMPUTES DAMPING FACTORS, EXTRAP ACTUALLY DOES THE DAMP.
C     DJ USED AS 3*L3 (3*L1*NORB, ACTUALLY)
C
      IF(ITER.EQ.0) GO TO 140
      IF(ITER.EQ.1) DEAVG = ZERO
      IF(ITER.EQ.2) DEAVG = ABS(DELE)
      IF(ITER.GE.3) DEAVG = (ABS(DELE)+ABS(DELE0)+PT2*DEAVG)/TWOPT2
C
      DMPTST = DENTOL
      IF ((DIFF0-DIFF) .LT. ZERO) DMPTST = DIFF0-DIFF
      IF(DAMPH  .AND.  ITER.GT.2  .AND.  ITDIIS.EQ.1)
     *    CALL DAMPD(DELE,DELE0,DEAVG,DAMP,DMPTST,DIFF,DIFFP,DMPTLC)
      IF (DAMP .LT. DMPCUT) DAMP = DMPCUT
C
      NDAF = 71
      IF((EXTRAH.OR.DAMPH)  .AND.  ITDIIS.EQ.1)
     *    CALL EXTRAP(DELE,DAMP,DAMP0,TRANS,DJ(1),DJ(LNOC+1),
     *                DJ(LNOC+LNOC+1),L1,LNOC,LNOC,NDAF,ITERV,1,2)
C
C     ----- ORTHONORMALIZE THE ORBITALS -----
C
      CALL ORTHO(Q,S,TRANS,SCR,L0,L0,L1,L2,L1)
      CALL TFSQB(TRANS,Q,SCR,L0,L1,L1)
      CALL DAWRIT(IDAF,IODA,TRANS,L3,15,0)
C
  140 CONTINUE
      IF(OUT) THEN
         WRITE(IW,*) 'CURRENT MO-S'
         CALL PREV(TRANS,E,L1,L1,L1)
      END IF
C
C     ----- SET UP THE DATA FOR THE VIRIAL CALCULATION -----
C     DJ USED AS L2
C     H  USED AS L2
C
      CALL VIRSET(TRANS,H,H2,DJ,L1,L2)
C
C     ----- FORM THE J AND K MATRICES -----
C     DJ USED AS NHAM*L2
C     H  USED AS NCORES*L2 + 2*(NHAM-NCORES)*L2
C
      CALL JKFORM(JKUNIT,TRANS,H,DJ,XP,XK,IXPK,NIJS,NKLS,IA,
     *            SCR,DSH,XINTS,GHONDO,DDIJ,L0,L1,L2,MXFOCK,
     *            NPASS,NINTMX,NOPK,INTTYP,NSH2,MAXG,SCHWRZ,
     *            DIRSCF,NINT,NSCHWZ,OUT)
      CALL TSECND(TIMDUM)
      TIMEF = TIMEF +TIMDUM - TIMITT
      TIMITT = TIMDUM
C
C     ----- DETERMINE THE OPTIMAL CI COEFFICIENTS -----
C     DJ USED AS 2*L2
C     H  USED AS 2*L2
C
      IF (.NOT.DIRSCF) DSKWRK = .TRUE.
      CALL OPTCI(JKUNIT,H,DJ,SCR,IWRK,IA,L1,L2,OUT)
      CALL TSECND(TIMDUM)
      TIMEG = TIMEG + TIMDUM - TIMITT
      TIMITT = TIMDUM
C
C     ----- EVALUATE THE FOCK MATRICES -----
C     DJ USED AS 2*L2 (DK IS SECOND PART OF DJ)
C     H  USED AS L2
C
      CALL GVBHAM(JKUNIT,DJ,DK,H,IA,L1,L2,OUT)
C
C     ----- SECOND ORDER GVB ORBITAL OPTIMIZATION -----
C
      SOSTART=SOSCF .AND. ITER.GE.NSO
      IF (SOSTART) THEN
         CALL GVBFCK(JKUNIT,H,DJ,S,TRANS,SCR,WRK3,WRK2,WRK1,DMPDII,
     *               DODIIS,L1,L2,L3,OUT,SOSCF)
         CALL GVBGRD(H,GRAD,L0,L2,NPR,ORBGRD)
         IF(ORBGRD.LT.SOGTOL .OR. ITSO.GT.0) THEN
            IF(ITSO.EQ.0) THEN
               EXTRAH=.FALSE.
               IF(MASWRK) WRITE(IW,9200)
            END IF
C
            ITSO = ITSO+1
            CALL GVBHES(JKUNIT,HESS,WRK1,DIAGH,DIAGJ,DIAGK,
     *                  L0,L2,NPR,NHAM)
            CALL SONEWT(HESS,GRAD,PGRAD,DISPL,WRK1,WRK2,WRK3,
     *                  WRK4,WRK5,WRK6,ORBGRD,NPR,ITSO,NFT15)
            CALL GVBTRN(DISPL,TRANS,WRK1,SCR,L0,L1,NPR,
     *                  SQCDF,RSTRCT)
            CALL DAWRIT(IDAF,IODA,WRK1,L3,61,0)
            CALL DCOPY(NPR,GRAD,1,PGRAD,1)
            GO TO 290
         END IF
      END IF
C
C     ----- SINGLE FOCK MATRIX ORBITAL OPTIMIZATION -----
C     FOCK OPERATOR IS RETURNED IN -H-
C     FOR DIIS, ERROR MATRIX IS FORMED AT -WRK3-
C            AND MO BASIS FOCK MATRIX IS IN -DJ-
C
      IF(NPAIR.LE.1) THEN
         CALL GVBFCK(JKUNIT,H,DJ,S,TRANS,SCR,WRK3,WRK2,WRK1,DMPDII,
     *               DODIIS,L1,L2,L3,OUT,SOSCF)
C
C        ----- POSSIBLE DIIS INTERPOLATION -----
C        IF DIIS TURNS ON, ITDIIS WILL BE .GT. 1
C        INTERPOLATED AO FOCK MATRIX MUST BE TRANSFORMED TO MO BASIS
C
         IF(DODIIS) THEN
            DSKWRK = SVDSKW
            CALL DIIS(GVB,IW,ITDIIS,Q,H,DMY,WRK3,DMY,WRK1,ADIIS,XDIIS,
     *                PDIIS,BDIIS,IODIIS,SCR,L1,L2,L3,MAXIT1,MAXIT2,
     *                MAXIO,ERDIIS,NOTOPN,MASWRK)
            CALL DCOPY(L2,H,1,DJ,1)
            CALL TFTRI(H,DJ,TRANS,SCR,L0,L1,L1)
         END IF
C
C        ----- LEVEL SHIFTING -----
C
         IF(VSHIFT .AND. ITER.GT.0  .AND.  ITDIIS.EQ.1)
     *         CALL LVLGVB(H,NCO,2,L0,DELE,DELE0,ITERV,1)
C
C        ----- DIAGONALIZE FOCK MATRIX -----
C
         IERR = 0
         CALL GLDIAG(L1,L0,L0,H,SCR,EIG,TRANS,IERR,IWRK)
         IF (IERR.NE.0) THEN
            IF(MASWRK) WRITE(IW,*) 'GVBITR,1: FAILED TO DIAGONALIZE F'
            CALL ABRT
         END IF
C
C        ----- SELECT OCCUPIED ORBITALS -----
C
         DO 150 IH = 1,NHAM
  150    CALL SELGV1(TRANS,EIG,L1,IH,RSTRCT)
         IF (OUT) WRITE(IW,*) 'ROTATION TO NEW ORBITALS'
         IF (OUT) CALL PREV(TRANS,EIG,L1,L1,L1)
         CALL DCOPY(L3,TRANS,1,DJ,1)
         CALL DCOPY(L1,EIG,1,E,1)
         CALL TSECND(TIMDUM)
         TIMEM = TIMEM + TIMDUM - TIMITT
         TIMITT = TIMDUM
      ELSE
C
C     ----- CONVENTIONAL GVB-ONE ORBITAL OPTIMIZATION -----
C     ----- DO THE 2X2 ROTATIONS AMONG OCCUPIED ORBITALS ----
C     DJ USED AS L3
C     H  USED AS 2*L2 (H2 IS SECOND PART OF H)
C
         ROTWT = ZERO
         IF (ROTDMP) ROTWT = DAMP
         ROTWT = ONE/(ONE+ROTWT)
         CALL ROTB(JKUNIT,DJ,H,H2,ROTWT,IA,L1,L2,OUT)
         CALL TSECND(TIMDUM)
         TIMEM = TIMEM + TIMDUM - TIMITT
         TIMITT = TIMDUM
C
C     ----- POSITION TO THE MO FOCK OPERATORS ----
C
         CALL SEQREW(JKUNIT)
         CALL SEQADV(JKUNIT)
         NHAM1 = NHAM
         IF(NCO.GT.0) NHAM1 = NHAM - 1
         IF(NHAM1.GT.0) THEN
            DO 180 I = 1,NHAM1
               CALL SEQADV(JKUNIT)
               CALL SEQADV(JKUNIT)
  180       CONTINUE
         END IF
C
C     ----- PERFORM THE OCBSE STEP TO MIX IN VIRTUALS -----
C
         DO 260 IHAM = 1,NHAM
            CALL SQREAD(JKUNIT,H,L2)
C
C     ----- ZERO OUT THE MATRIX ELEMENTS CONNECTING THE ORBITALS OF
C           THE IHAM'TH FOCK OPERATOR TO THE ORBITALS OF THE OTHER
C           FOCK OPERATORS, THE LOOP SHIFTS DIAGONAL ELEMENTS OF THE
C           OTHER FOCK OPERATORS OUT OF THE INTERACTION RANGE -----
C
            DO 220 I = 1,L1
               EIG(I) = 1000-3*I
  220       CONTINUE
            CALL ZEROH(H,EIG,IHAM,IA,L1,L2)
C
C     ----- INVOKE LEVEL SHIFTING IF NEEDED -----
C
            IF(NCO.GT.0  .AND.  IHAM.EQ.1) THEN
               NSHFT = NORB
               NCALL = 1
            ELSE
               NSHFT = 0
               NCALL = 0
            END IF
            IF(VSHIFT .AND. ITER.GT.0)
     *         CALL LVLGVB(H,NSHFT,0,L0,DELE,DELE0,ITERV,NCALL)
C
C     ----- DIAGONALIZE THE CURRENT FOCK MATRIX -----
C
            IERR = 0
            CALL GLDIAG(L1,L1,L1,H,SCR,EIG,TRANS,IERR,IWRK)
            IF (IERR.NE.0) THEN
               IF(MASWRK) WRITE(IW,*)
     *             'GVBITR,2: FAILED TO DIAGONALIZE F'
               CALL ABRT
            END IF
C
C     ----- SELECT THE ORBITALS WHICH LOOK MOST LIKE CURRENT ITER ----
C     ----- UPDATE THE ORBITALS WITH THE NEW OCBSE CONTRIBUTIONS -----
C
            CALL SELGV1(TRANS,EIG,L1,IHAM,RSTRCT)
            CALL UPOCBS(DJ,TRANS,E,EIG,L1,IHAM)
  260    CONTINUE
         CALL TSECND(TIMDUM)
         TIMEO = TIMEO + TIMDUM - TIMITT
         TIMITT = TIMDUM
      END IF
C
C     ----- SCHMIDT ORTHONORMALIZE THE ORBITAL ROTATION MATRIX -----
C     ----- CALCULATE EXTENT OF ORBITAL ROTATION ----
C
      DSKWRK = SVDSKW
      CALL SCHORT(DJ,L0,L0,L1)
      CALL GVBCVG(DJ,E,L1,L0,SQCDF,OUT)
C
C     ----- SAVE THE ROTATION MATRIX -----
C     ----- TRANSFORM THE ORBITALS BACK TO THE AO BASIS -----
C     ----- CLEAN UP THE ORBITALS -----
C
      CALL DAWRIT(IDAF,IODA,DJ,L3,61,0)
      CALL DAREAD(IDAF,IODA,TRANS,L3,15,0)
      CALL TFSQB(DJ,TRANS,SCR,L0,L1,L1)
      CALL DCOPY(L3,DJ,1,TRANS,1)
C
C           AT THIS POINT, ALL CONVERGENCE METHODS SHOULD HAVE
C           PLACED THE ORBITALS IN THE AO BASIS IN -TRANS-
C
  290 CONTINUE
      CALL CLENMO(TRANS,L1,L0,TOLZ,TOLE,IW,OUT)
C
C     ----- CALCULATE THE REDUCED DENSITY MATRICES -----
C
      CALL REDDEN(TRANS,DJ,DK,L1,L2)
C
C     ----- COMBINE THE ALPHA AND BETA INTO ONE DENSITY MATRIX -----
C     GET THE DENSITY DIFFERENCE AND SAVE THE NEW DENSITY MATRIX
C
      DO 300 I = 1,L2
         DJ(I) = DJ(I) + DK(I)
  300 CONTINUE
      IF(ITER.GT.0) THEN
         CALL DAREAD(IDAF,IODA,DK,L2,62,0)
      ELSE
         CALL VCLR(DK,1,L2)
      END IF
      DIFFP = DIFF
      CALL DDIFF(DK,DJ,L2,DIFF)
      CALL DAWRIT(IDAF,IODA,DJ,L2,62,0)
C
C        ----- SOLVENT ENERGY CORRECTIONS -----
C     NOTE: ONSAGER MODEL PUT FREE ENERGY CORRECTIONS INTO -EN- ALREADY
C
      EHF = EHF - EEFP
C
      IF(IPCFP.OR.IPCM.EQ.1) THEN
          IF(IPCMIT.EQ.0) THEN
             EHF=EHF-HALF*(PB+PX-PC)
          ELSE
             EHF=EHF+PX-HALF*ENPCM
          END IF
      END IF
      IF(IPCM.EQ.1) THEN
         EHF=EHF+(GCAVP+GDISP+GREP)/627.509541D+00
      END IF
C
C     SVP MODEL'S FREE ENERGY CORRECTION IS
C     +HALF*ATMSVP SINCE IT HAS NOT YET BEEN ACCOUNTED FOR, BUT
C     -HALF*ELCSVP SINCE EHF ALREADY IMPLICITLY CONTAINS ELCSVP
C
      IF(ISVP.EQ.1) THEN
         EHF=EHF+HALF*(ATMSVP-ELCSVP)
      END IF
C
C     ----- PRINT ITERATION RESULTS -----
C
      EHF   = EHF + EHNCO
      EHF0  = EHF
      DIFF0 = DIFF
      ETOT0 = ETOT
      ETOT  = EHF+EN
      DELE0 = DELE
      DELE  = ETOT - ETOT0
      EK  = EK+EK
      VIR = (ETOT-EK)/(TWO*ETOT)
      IF (MASWRK) THEN
         IF(SOSCF) THEN
            ERRORC = ORBGRD
         ELSE
            ERRORC = ERDIIS
         END IF
      IF(DAMPH  .OR.  VSHIFT) THEN
         IF(DIRSCF) THEN
            WRITE(IW,9090)
     *      ITER,ICALP,ETOT,DELE,SQCDF,ERRORC,RRSHFT,DAMP,NINT,NSCHWZ
         ELSE
            WRITE(IW,9090)
     *      ITER,ICALP,ETOT,DELE,SQCDF,ERRORC,RRSHFT,DAMP
         END IF
      ELSE
         IF(DIRSCF) THEN
            WRITE(IW,9100)
     *      ITER,ICALP,ETOT,DELE,SQCDF,ERRORC,NINT,NSCHWZ
         ELSE
            WRITE(IW,9100)
     *      ITER,ICALP,ETOT,DELE,SQCDF,ERRORC
         END IF
      END IF
      END IF
      CALL FLSHBF(IW)
      ICALP = ICALP+1
C
C     ----- CHECK FOR CONVERGENCE -----
C
      CVDENS = (DIFF.LT.DENTOL  .AND.  ABS(DELE).LT.ETOL2)  .OR.
     *         (DIFF.LT.DENTOL/FIVE)
      CVENGY = ABS(DELE).LT.ENGTOL  .AND.  DIFF.LT.DTOL2
      IF(DODIIS) CVDIIS = ERDIIS.LT.DIITOL  .AND.  DIFF.LT.DTOL2
      CVGED  = CVGING  .AND.  (CVDENS.OR.CVENGY.OR.CVDIIS)
      CVGING = CVDENS.OR.CVENGY.OR.CVDIIS
C
      IF(ISVP.EQ.1) THEN
         CVGED  = CVGED   .AND.  CVGSVP
         CVGING = CVGING  .AND.  CVGSVP
         CVGSVPNR = DIFF  .LT.  100.0D+00 * DENTOL
         CVGSVPDN = CVGING
         CVGSVPDN1 = CVGING
      END IF
C
      IF(CVGED) THEN
         IF(CVDENS) THEN
            IF (MASWRK) WRITE(IW,9124)
            GO TO 360
         END IF
         IF(CVDIIS) THEN
            IF (MASWRK) WRITE(IW,9126)
            GO TO 360
         END IF
         IF(CVENGY) THEN
            IF (MASWRK) WRITE(IW,9128)
            GO TO 360
         END IF
      END IF
      ITER = ITER + 1
C
C     ----- CHECK FOR TIME REMAINING -----
C
      CALL TSECND(TIME1)
      TLEFT = TIMLIM - TIME1
      TIME0 = TIME1
      IF (TLEFT .LT. TWO*(TLEFTS-TLEFT)/ITER) THEN
         IF (MASWRK) WRITE (IW,9088)
         ETOT = ZERO
         EHF = -EN
         GO TO 360
      END IF
C
      IF (ITER .LE. MAXIT) GO TO 100
C
C           *********************
C     ***** END OF GVB ITERATIONS *****
C           *********************
C
      IF (MASWRK) WRITE (IW,9108)
      ETOT = ZERO
      EHF = -EN
  360 CONTINUE
C
C     ----- CANONICALIZE ORBITALS ----
C
      IF (SOSCF) THEN
C               DIAGONALIZE COMBINED FOCK MATRIX
         CALL GVBDIA(H,TRANS,E,SCR,IWRK,L0,L1,L2,L1)
C               POSSIBLY RESTRICT ORBITAL INTERCHANGES
         DO 370 IHAM = 1,NHAM
  370    CALL SELGV1(TRANS,E,L1,IHAM,RSTRCT)
C               GET TOTAL ORBITAL ROTATION (FINAL SOSCF X CANON)
         CALL DAREAD(IDAF,IODA,DJ,L3,61,0)
         CALL MRARBR(DJ,L0,L0,L0,TRANS,L1,L0,DK,L1)
         CALL DCOPY(L1,DK,1,TRANS,1)
         CALL DAWRIT(IDAF,IODA,TRANS,L3,61,0)
C               BACK TRANSFORM TO GET ORBITALS IN AO BASIS
         CALL DAREAD(IDAF,IODA,DJ,L3,15,0)
         CALL TFSQB(TRANS,DJ,SCR,L0,L1,L1)
         CALL DAWRIT(IDAF,IODA,TRANS,L3,15,0)
C               GENERATE AND SAVE DENSITY FOR THESE ORBS
         CALL REDDEN(TRANS,DJ,DK,L1,L2)
      END IF
C
C     ----- PRINT FINAL RESULTS -----
C
      IF (MASWRK) WRITE(IW,9130) ETOT,ITER
      IF(MPCTYP.NE.NONE  .AND.  ETOT.NE.ZERO) THEN
         DELTAH = ETOT*627.52791D+00 + ATHEAT
         IF(MASWRK) WRITE(IW,9132) DELTAH
      END IF
C
C     ----- CALCULATE THE FIVE TIME STEPS PER ITERATION -----
C
      TIMEF = TIMEF/ ITER
      TIMEG = TIMEG/ ITER
      TIMEM = TIMEM/ ITER
      TIMEO = TIMEO/ ITER
      IF(SOME) WRITE (IW,9148) NPASS,TIMEF,TIMEG,TIMEM,TIMEO
C
C     ----- ADJUST THE ENERGIES OF THE PAIR ORBITALS ONLY -----
C     I'M NOT SURE THAT MULTIPLYING BY F(I) MAKES SENSE, BUT
C     THAT'S THE WAY THE PROGRAM WAS (EXCEPT IT ALSO MULTIPLIED
C     THE ENERGIES BY 2) - MWS.
C
      IF (NPAIR.GT.0) THEN
         NSKIP = 0
         IF(NCO.GT.0) NSKIP = 1
         NSKIP = NSKIP + NSETO
         DO 400 I = 1,NORB
            IF(NCONF(I).LE.NSKIP) GO TO 400
            E(I) = E(I) * F(NCONF(I))
  400    CONTINUE
      END IF
C
C     ----- FINISH THE TCSCF HAMILTONIAN MATRIX -----
C
      IF(NPAIR.EQ.1) CALL GVBTCH(ETOT)
C
C     ----- SAVE THE FINAL EIGENVALUES AND VECTORS -----
C     ALPHA AND BETA SETS SAVED FOR COMPATABILITY WITH PROPERTIES CODE
C
      CALL DAWRIT(IDAF,IODA,TRANS,L3,15,0)
      CALL DAWRIT(IDAF,IODA,TRANS,L3,19,0)
      CALL DAWRIT(IDAF,IODA,E,L1,17,0)
      CALL DAWRIT(IDAF,IODA,E,L1,21,0)
C
C     ----- SAVE THE LAGRANGIAN MULTIPLIERS ON THE DAF -----
C     DJ IS USED AS NORB*NORB
C     H  IS USED AS L2
C
      CALL DAREAD(IDAF,IODA,TRANS,L3,61,0)
      CALL GVBLAG(JKUNIT,L0,L1,NORB,L2,H,DJ,TRANS,IA)
      L2NORB = NORB*NORB
      CALL DAWRIT(IDAF,IODA,DJ,L2NORB,402,0)
C
C         SOLVATION OUTPUT FOR PCM MODEL
C
      IF(IPCFP) THEN
         IPCM=1
         IEFP=1
         IPCFP=.FALSE.
      END IF
C
      QTT=0.0D+00
      DO I=0,NTS-1
C        QT=QSE(I)+QSN(I)+Q_FS(I)+Q_IND(I)
         QT=X(LQSE+I)+X(LQSN+I)+X(LQFS+I)+X(LQIND+I)
         QTT=QTT+QT
      ENDDO
      IF(IPCM.EQ.1.AND.IPCMIT.EQ.1) THEN
        IF(MASWRK) WRITE(IW,*) ' '
        IF(MASWRK) WRITE(IW,'(A10,F18.10)')'PCM ASC=',QTT
      END IF
      IF(IPCM.EQ.1) THEN
        IF(MASWRK) WRITE(IW,'(A10,F18.10)')
     *      'GCAVP  =',GCAVP/627.509541D+00
        IF(MASWRK) WRITE(IW,'(A10,F18.10)')
     *      'GDISP  =',GDISP/627.509541D+00
        IF(MASWRK) WRITE(IW,'(A10,F18.10)')
     *      'GREP   =',GREP/627.509541D+00
      END IF
C
      IF(MASWRK  .AND.  IPCM.EQ.1  .AND.  IPCMIT.EQ.0) THEN
         WRITE(IW,*)
         WRITE(IW,*) '-- INDUCED PCM ELECTRONIC CHARGES --'
         IF(ICOMP.EQ.0) WRITE(IW,1010) QET, TCH
         IF(ICOMP.EQ.2) WRITE(IW,1100) QET, TCH, QETN
         IF(ICOMP.EQ.3) WRITE(IW,1110) QET, TCH, QESC, QETN
         WRITE(IW,*) '--------------------------------------'
      END IF
C
C         SOLVATION OUTPUT FOR SVP MODEL
C
      IF(ISVP.EQ.1  .AND.  MASWRK) THEN
         CALL SVPOUT(ETOT,ECHIP,NPRINT,'GVBITR')
      END IF
C
C         SOLVATION OUTPUT FOR ONSAGER OR EFP MODELS
C
      IF ((IZRF.EQ.1 .OR. IEFP.EQ.1) .AND.  MASWRK)
     *     WRITE(IW,9135) TMUX,TMUXD,TMUY,TMUYD,TMUZ,TMUZD
C
C     ----- MODIFIED VIRTUAL ORBITALS -----
C     NOTE THAT D AND FMO ARE EQUIVALENT BY CALL (BOTH IN H2)
C
      IF(MVOQ.GT.0  .AND.  EXETYP.NE.CHECK) THEN
         NFROZ = NORB
         CALL MVOS(TRANS,E,H,H2,SCR,IA,DJ,H2,XP,IXPK,
     *             DSH,GHONDO,DDIJ,XINTS,DIRSCF,SCHWRZ,
     *             NOPK,NINTMX,INTTYP,L0,L1,L2,L3,NSH2,MAXG,
     *             NFROZ,MVOQ,NE)
      END IF
C
C     --- GENERAL CORRELATION ADAPTED VIRTUAL ORBITALS ---
C     STORE SYMMETRIC -DCORE-, -DVAL- AT -WRK1-, WRK2-
C     STORE SYMMETRIC -JCORE-, -JVAL- AT -H-, -H2-
C     STORE SYMMETRIC -KCORE-, -KVAL- AT -WRK3-, -S-
C
      IF(IACAVO.EQ.1  .AND.  EXETYP.NE.CHECK) THEN
         NFROZ = NORB
         SVDSKW = DSKWRK
         DSKWRK = .TRUE.
         CALL ACAVO(PACAVO,TRANS,E,WRK1,WRK2,
     *              H,H2,WRK3,S,
     *              Q,SCR,IWRK,XP,IXPK,
     *              NOPK,NINTMX,L0,L1,L2,L3,NFROZ)
         DSKWRK = SVDSKW
      END IF
C
C     ----- PRINT AND PUNCH THE ORBITALS -----
C
      CALL DAREAD(IDAF,IODA,TRANS,L3,15,0)
      CALL GVBOUT(TRANS,S,Q,IWRK,SCR(1,1),SCR(1,3),E,SCR(1,4),
     *            L0,L1,L2,L3)
C
  700 CONTINUE
      IF (.NOT.DIRSCF) DSKWRK = .TRUE.
      CALL SEQCLO(JKUNIT,'KEEP')
      IF(.NOT.NOTOPN) CALL RACLOS(20,'DELETE')
      NOTOPN = .TRUE.
      DSKWRK = SVDSKW
C
C        RETURN SOLVENT MODEL MEMORY
C
      IF (IZRFLS .OR. IEFP.EQ.1 .OR.
     *    IPCM.EQ.1 .OR. NBUFMO.GT.0  .OR.  ISVP.EQ.1) THEN
          CALL RETFM(NEED)
      END IF
C
      RETURN
C
 9010 FORMAT(/5X,'RSTRCT FORCED TRUE BY THE PRESENCE OF DEGENERATE',
     *           ' OPEN SHELLS.')
 9017 FORMAT(5X,'SOSCF WILL OPTIMIZE',I8,' ORBITAL ROTATIONS, SOGTOL=',
     *          F8.3)
 9018 FORMAT(/5X,'MAXIT=',I4,'   NPUNCH=',I2,'   SQCDF TOL=',1P,E10.4/
     *       5X,'NUCLEAR ENERGY=',0P,F20.10/5X,'EXTRAP=',L1,
     *     '   DAMP=',L1,'   SHIFT=',L1,'   RSTRCT=',L1,'   DIIS=',L1,
     *     '  SOSCF=',L1)
 9020 FORMAT(1X,'SOLVATION MODEL FOR GVB REQUIRES',I10,
     *    ' ADDITIONAL WORDS OF MEMORY')
 9030 FORMAT(/1X,'DIRECT SCF CALCULATION, SCHWRZ=',L1,'   FDIFF=',L1)
 9040 FORMAT(/' ITER EX     TOTAL ENERGY       E CHANGE        ',
     *        'SQCDF       ORB. GRAD       VIR. SHIFT       DAMPING')
 9045 FORMAT(/' ITER EX     TOTAL ENERGY       E CHANGE        ',
     *        'SQCDF       DIIS ERROR      VIR. SHIFT       DAMPING')
 9050 FORMAT(/107X,'NONZERO    BLOCKS'/
     *        ' ITER EX     TOTAL ENERGY       E CHANGE        ',
     *        'SQCDF       ORB. GRAD       VIR. SHIFT       DAMPING',
     *        '        INTEGRALS   SKIPPED')
 9055 FORMAT(/107X,'NONZERO    BLOCKS'/
     *        ' ITER EX     TOTAL ENERGY       E CHANGE        ',
     *        'SQCDF       DIIS ERROR      VIR. SHIFT       DAMPING',
     *        '        INTEGRALS   SKIPPED')
 9060 FORMAT(/' ITER EX     TOTAL ENERGY       E CHANGE        ',
     *        'SQCDF       ORB. GRAD ')
 9065 FORMAT(/' ITER EX     TOTAL ENERGY       E CHANGE        ',
     *        'SQCDF       DIIS ERROR')
 9070 FORMAT(/75X,'NONZERO    BLOCKS'/
     *        ' ITER EX     TOTAL ENERGY       E CHANGE        ',
     *        'SQCDF       ORB. GRAD       INTEGRALS   SKIPPED')
 9075 FORMAT(/75X,'NONZERO    BLOCKS'/
     *        ' ITER EX     TOTAL ENERGY       E CHANGE        ',
     *        'SQCDF       DIIS ERROR      INTEGRALS   SKIPPED')
 9090 FORMAT(1X,2I3,F19.9,F16.9,2F14.9,2F16.9,I15,I10)
 9100 FORMAT(1X,2I3,F19.9,F16.9,2F14.9,I15,I10)
 9088 FORMAT(/1X,'SCF IS UNCONVERGED, TOO LITTLE TIME LEFT')
 9108 FORMAT(/1X,'SCF IS UNCONVERGED, TOO MANY ITERATIONS')
 9124 FORMAT(/10X,17(1H-)/10X,17HDENSITY CONVERGED/10X,17(1H-))
 9126 FORMAT(/10X,14(1H-)/10X,14HDIIS CONVERGED/10X,14(1H-))
 9128 FORMAT(/10X,16(1H-)/10X,16HENERGY CONVERGED/10X,16(1H-))
 9130 FORMAT(/1X,'FINAL GVB ENERGY IS',F20.10,' AFTER',I4,' ITERATIONS')
 9132 FORMAT(/1X,'HEAT OF FORMATION IS',F15.5,' KCAL/MOL')
 9135 FORMAT(/10X,'INDUCED DIPOLE',5X,'ATOMIC UNITS',10X,'DEBYE',
     *       /10X,'      X(IND)  ',5X,F12.5,5X,F12.5,
     *       /10X,'      Y(IND)  ',5X,F12.5,5X,F12.5,
     *       /10X,'      Z(IND)  ',5X,F12.5,5X,F12.5)
 9148 FORMAT(/,10X,28(1H-),/,10X,28HSCF STATISTICS PER ITERATION,/,10X,
     *     28(1H-),/,10X, 28HNUMBER OF INTEGRAL PASSES   ,I5,/,10X,
     *     23HFOCK FORMATION TIME    ,F10.3,/,10X,
     *     23HGEMINAL OPT TIME       ,F10.3,/,10X,
     *     23HMIXORB  OPT TIME       ,F10.3,/,10X,
     *     23HOCBSE   OPT TIME       ,F10.3)
 9200 FORMAT(10X,15(1H-),'START SECOND ORDER SCF',15(1H-))
C
 1010 FORMAT(1X,'CALCULATED ',F10.5,' THEORETICAL ',
     * F10.5,' NOT RENORMALIZED',F10.5)
 1100 FORMAT(1X,'CALCULATED ',F10.5,' THEORETICAL ',
     * F10.5,' RENORMALIZED',F10.5)
 1110 FORMAT(1X,'CALCULATED',F10.5,' (THEOR= ',
     * F10.5,')  ESCAPED',F10.5,' FINAL', F10.5)
      END
C*MODULE GVB     *DECK GVBLAG
      SUBROUTINE GVBLAG(JKUNIT,L0,L1,NORBX,L2,H,DJ,TRANS,IA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION H(L2),DJ(NORBX,*),TRANS(L1,*),IA(L1)
C
      COMMON /GVBENG/ EK,VIR,TKIN(25),TPLUSV(25),ASYM
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- CALCULATE THE LAGRANGIAN MULTIPLIER MATRIX -----
C
C     SKIP PAST THE J AND K MATRICES TO THE FOCK MATRICES
C
      CALL SEQREW(JKUNIT)
      CALL SEQADV(JKUNIT)
      NHAM1=NHAM
      IF(NCO.GT.0) NHAM1=NHAM-1
      DO 100 I = 1,NHAM1
         CALL SEQADV(JKUNIT)
         CALL SEQADV(JKUNIT)
  100 CONTINUE
C
C     ----- CONSTRUCT THE LAGRANGIAN MATRIX IN MO BASIS -----
C     UPDATE THE FOCK OPERATORS IN THE MO BASIS, FOR THE
C     FINAL SMALL ROTATION MATRIX MADE BY THE LAST ITERATION.
C
C     A CASE WHERE DEGENERATE ORBITALS OCCUR IN DIFFERENT
C     SHELLS PROBABLY HAS A LARGE ANGLE OF ROTATION BETWEEN
C     THEM IN THE LAST ITERATION.  ROTATING THE FOCK MATRICES
C     DOES NOT GIVE THE CORRECT LAGRANGIAN IN THIS CASE, AT
C     LEAST FOR CASES LIKE THE 1-DELTA STATES (BOTH XY AND
C     XX-YY COMPONENTS) ARISING FROM PI**2 AND COMPUTED WITH
C     NOSYM=1 SO THAT THE ENERGY IS CORRECT.  SINCE WE DO NO
C     KNOW HOW TO FIX IT, WE GENERATE A WARNING MESSAGE BELOW,
C     AND ALSO PLACE A NOTE IN THE MANUAL.  NOTE THAT SOSCF
C     CONVERGENCE DOES NOT HAVE A ROTATION BETWEEN DEGENERATE
C     OPEN SHELLS, AND SO GIVES A GOOD LAGRANGIAN BELOW
C
      CALL VCLR(DJ,1,NORB*NORB)
C
      DO 290 IHAM = 1,NHAM
C
C           READ THIS SHELL'S FOCK MATRIX IN MO BASIS
C
         CALL SQREAD(JKUNIT,H,L2)
         IF (GOPARR.AND.(.NOT.DSKWRK))
     *       CALL DDI_BCAST(1025,'F',H,L2,MASTER)
C
C            ROTATE IT, AND SET EACH COLUMN OF THE LAGRANGIAN
C            FOR ANY ORBITAL OCCUPIED IN THIS SHELL.
C
         DO 280 I = 1,NORB
            IF (NCONF(I).NE.IHAM) GO TO 280
            DO 250 J = 1,NORB
               EPS = ZERO
               DO 240 K = 1,L0
                  DO 230 L = 1,L0
                     KL = IA(K)+L
                     IF (L .GT. K) KL = IA(L)+K
                     EPS = EPS + TRANS(K,J)*H(KL)*TRANS(L,I)
  230             CONTINUE
  240          CONTINUE
               DJ(J,I) = EPS
  250       CONTINUE
  280    CONTINUE
  290 CONTINUE
C
C     ----- CONVERT TO THE GENERALIZED FOCK OPERATOR -----
C
      DO 360 I = 1,NORB
         NX = NCONF(I)
         DO 320 J = 1,NORB
            DJ(J,I) = DJ(J,I)*F(NX)
  320    CONTINUE
         DO 340 J = 1,NORB
            DJ(J,I) = DJ(J,I)+DJ(J,I)
  340    CONTINUE
  360 CONTINUE
C
C     ----- FIND THE LARGEST ASYMMETRY -----
C
      ASYM = ZERO
      DO 420 I = 1,NORB
         DO 410 J = 1,I
            DIFF = ABS(DJ(I,J)-DJ(J,I))
            IF (DIFF .GT. ASYM) ASYM = DIFF
  410    CONTINUE
  420 CONTINUE
C
C     ----- WARN OF POTENTIAL PROBLEMS ----
C
      IF(ASYM.GT.1.0D-05  .AND.  MASWRK) WRITE(IW,9000)
      RETURN
C
 9000 FORMAT(//1X,58(1H*)/
     *  1X,'WARNING * * * THE LAGRANGIAN MATRIX IS NOT VERY SYMMETRIC.'/
     *  1X,'THIS MAY BE DUE TO A FAILURE TO CONVERGE THE GVB-SCF.'/
     *  1X,'IT MIGHT BE DUE TO USING DIIS WITH DEGENERATE OPEN SHELLS.'/
     *  1X,'IN EITHER CASE, GRADIENTS AND HESSIANS WILL BE --WRONG--.'/
     *  1X,'REMEDIES: TRY TO CONVERGE THE SCF, OR USE SOSCF METHOD.'/
     *  1X,58(1H*)//)
      END
C*MODULE GVB     *DECK GVBOUT
      SUBROUTINE GVBOUT(TRANS,S,Q,IWRK,GI,INORDR,E,WRK,L0,L1,L2,L3)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      LOGICAL SOME,LOTS,GOPARR,DSKWRK,MASWRK
C
      CHARACTER*8 GWRD,EWRD
C
      DIMENSION TRANS(L1,L1),S(L2),Q(L3),IWRK(L1),GI(L1,2),INORDR(L1),
     *          E(L1),WRK(L1)
      DIMENSION TIMSTR(3)
C
      COMMON /CONV  / DENTOL,EN,ETOT,EHF,EHFO,DIFF,ITER,ICALP,ICBET
      COMMON /GVBCI / CIHAM(91),CILOW(12),JPAIR,KONE,KTWO,KCORB(2,12)
      COMMON /GVBENG/ EK,VIR,TKIN(25),TPLUSV(25),ASYM
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /MEXOPT/ MEXSKPGES,MEXSTATE
      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 /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
C
      DATA AMEX/8HMEX     /
C
C     ----- ROHF-GVB OUTPUT ROUTINE -----
C
      SOME = NPRINT.NE.-5 .AND. MASWRK
      LOTS = NPRINT.EQ.8 .AND. MASWRK
C
C     ----- PRINT PAIR INFORMATION -----
C
      IF(NPAIR.GT.0) THEN
         IF (MASWRK) THEN
            WRITE (IW,9108)
            WRITE (IW,9128)
         END IF
         DO 120 I = 1,NPAIR
            SAB = (CICOEF(1,I)+CICOEF(2,I))/(CICOEF(1,I)-CICOEF(2,I))
            SAB = ABS(SAB)
            OCC1 = TWO * CICOEF(1,I)**2
            OCC2 = TWO * CICOEF(2,I)**2
            IF (MASWRK) WRITE (IW,9148)
     *         I,(KCORB(J,I),J=1,2),(CICOEF(J,I),J=1,2),
     *         OCC1,OCC2,SAB,CILOW(I)
  120    CONTINUE
      END IF
C
C     ----- ORDER VIRTUAL ORBITALS BY ORBITAL ENERGY -----
C
      NVIR = L0-NORB
      IVIR = NORB+1
      CALL ORDERV(TRANS(1,IVIR),E(IVIR),INORDR,GI,NVIR,NVIR,L1)
C
      IF(MASWRK) WRITE (IW,9288) ASYM
C
C     ----- ENERGY COMPONENTS, AND EIGENVALUE PRINTOUT -----
C
      IF(LOTS) THEN
         WRITE (IW,9168) EK,EHF,EN,ETOT,VIR
         WRITE (IW,9048)
         DO 140 I = 1,NORB
            WRITE (IW,9068) I,E(I)
  140    CONTINUE
      END IF
C
C     ----- AND PRINT THE ORBITALS -----
C     NOTE THAT THE CALL TO SYMMOS WILL DESTROY THE -Q- MATRIX
C
      CALL SYMMOS(IWRK,Q,S,TRANS,WRK,L0,L1,L0,L1)
      CALL DAWRIT(IDAF,IODA,IWRK,L1,255,1)
C
      IF(SOME) THEN
         WRITE (IW,9088)
         CALL PREVS(TRANS,E,IWRK,L0,L1,L1,1)
      END IF
C
C     ----- PUNCH THE ORBITALS -----
C
      IF(MASWRK  .AND.  NPUNCH.GT.0) THEN
         CALL TMDATE(TIMSTR)
         WRITE(IP,8000) TIMSTR,TITLE,ETOT,EN,ITER
         DO 160 I=1,NPAIR
            INDEX = 2*I - 1
            GWRD = '        '
            EWRD = '        '
            IF(I.EQ.1)     GWRD=' $SCF   '
            IF(I.EQ.NPAIR) EWRD=' $END   '
            WRITE(IP,8020) GWRD,INDEX,CICOEF(1,I),CICOEF(2,I),EWRD
  160    CONTINUE
         IF(RUNTYP.NE.AMEX) THEN
            WRITE(IP,8010) ' $VEC   '
         ELSE
            IF(MEXSTATE.EQ.1) WRITE(IP,8010) ' $VEC1  '
            IF(MEXSTATE.EQ.2) WRITE(IP,8010) ' $VEC2  '
         END IF
         LPUN=L0
         IF(NPUNCH.EQ.1) LPUN=NORB
         CALL PUSQL(TRANS,LPUN,L1,L1)
         WRITE(IP,8010) ' $END   '
      END IF
C
C     ----- COMPUTE AND PRINT THE GVB PAIR ORBITALS -----
C     THESE ARE THE NONORTHOGONAL ORBITALS GODDARD LIKES TO PLOT.
C
      IF(SOME  .AND.  NPAIR .GT. 0) THEN
         WRITE (IW,9188)
         DO 360 I = 1,NPAIR
            C1 = CICOEF(1,I)
            C2 = CICOEF(2,I)
            M1 = KCORB(1,I)
            M2 = KCORB(2,I)
            ALPH = SQRT(-C1/C2)
            BET = ONE/SQRT(ONE+ALPH**2)
            ALPH = ALPH*BET
            RMAX1 = ZERO
            RMAX2 = ZERO
            IPH1 = 1
            IPH2 = 1
            DO 320 K = 1,L1
               GI(K,1) = ALPH*TRANS(K,M1) + BET*TRANS(K,M2)
               GI(K,2) = -ALPH*TRANS(K,M1)+BET*TRANS(K,M2)
               IF ( ABS(GI(K,1)) .GT. RMAX1) THEN
                  IPH1 = 1
                  RMAX1 = ABS(GI(K,1))
                  IF (GI(K,1) .LT. ZERO) IPH1 = -1
               END IF
               IF (ABS(GI(K,2)) .GT. RMAX2) THEN
                  IPH2 = 1
                  RMAX2 = ABS(GI(K,2))
                  IF (GI(K,2) .LT. ZERO) IPH2 = -1
               END IF
  320       CONTINUE
            DO 340 K = 1,L1
               IF (IPH1 .GT. 0) GO TO 340
               GI(K,1) = -GI(K,1)
               IF (IPH2 .GT. 0) GO TO 340
               GI(K,2) = -GI(K,2)
  340       CONTINUE
            WRITE (IW,9208) I
            CALL PRSQL(GI,2,L1,L1)
  360    CONTINUE
      END IF
C
C     ----- PRINT OUT THE FINAL F COEFFICIENT MATRIX -----
C
      IF(LOTS) THEN
         WRITE (IW,9228)
         DO 400 I = 1,NHAM
            WRITE (IW,9248) I,F(I)
  400    CONTINUE
      END IF
C
C     ----- PRINT THE LAGRANGIAN MATRIX -----
C     ALL PROCESS MUST DAREAD, BUT ONLY MASTER CAN PRINT
C
      IF(NPRINT.EQ.8) THEN
         L2NORB = NORB*NORB
         CALL DAREAD(IDAF,IODA,TRANS,L2NORB,402,0)
      END IF
      IF(LOTS) THEN
         WRITE (IW,9268)
         CALL PRSQ(TRANS,NORB,NORB,NORB)
      END IF
      RETURN
C
 8000 FORMAT('--- GVB ORBITALS --- GENERATED AT ',3A8/10A8/
     *       'E(GVB)=',F20.10,', E(NUC)=',F16.10,',',I5,' ITERS')
 8010 FORMAT(A8)
 8020 FORMAT(A8,' CICOEF(',I2,')=',F12.8,',',F12.8,A8)
C
 9048 FORMAT(/,10X,11(1H-),/,10X,11HEIGENVALUES,/,10X,11(1H-),/)
 9068 FORMAT(1X,I5,F20.12,4X,F20.12)
 9088 FORMAT(/10X,12(1H-)/10X,'EIGENVECTORS'/10X,12(1H-))
 9108 FORMAT(/10X,16(1H-)/10X,'PAIR INFORMATION'/10X,16(1H-))
 9128 FORMAT(1X,5X,'ORBITAL',3X,'CI COEFFICIENTS',4X,
     *          'OCCUPATION NUMBERS',3X,'GVB',6X,'ENERGY'/
     *       1X,'PAIR 1   2',5X,'ORB 1',5X,'ORB 2',6X,'ORB 1',
     *          5X,'ORB 2',4X,'OVERLAP  LOWERING')
 9148 FORMAT(1X,I2,2I4,2X,2F10.6,4F10.5)
 9168 FORMAT(/10X,
     +     20HKINETIC ENERGY       ,F16.10/ 10X,20HELECTRONIC ENERGY
     +     ,F16.10/ 10X,14HNUCLEAR ENERGY,6X,F16.10,5X,
     +     18HVIRIAL RATIO(V/2E),/, 10X,20HTOTAL  ENERGY        ,F16.10,
     +     8X,F10.7)
 9188 FORMAT(/,10X,11(1H-)/10X,11HGI ORBITALS/10X,11(1H-))
 9208 FORMAT(//,20X,5HPAIR ,I3)
 9228 FORMAT(/,10X,20(1H-),/,10X,20HF COEFFICIENT MATRIX,/,10X, 20(1H-),
     +     / )
 9248 FORMAT(3X,5X,I5,F20.12)
 9268 FORMAT(/,10X,22(1H-),/,10X,22HLAGRANGIAN MULTIPLIERS,/,10X,22(1H-)
     +     , /)
 9288 FORMAT(/5X,'THE MAXIMUM LAGRANGIAN ASYMMETRY IS ',1P,E13.7)
      END
C*MODULE GVB     *DECK GVBSET
      SUBROUTINE GVBSET(OUT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL OUT
C
      COMMON /GVBCI / CIHAM(91),CILOW(12),JPAIR,KONE,KTWO,KCORB(2,12)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
C
C     GENERATE SOME INITIAL PARAMETERS FOR THE GVB RUN
C
C     ----- GET TOTAL NUMBER OF FOCK OPERATORS -----
C
      NCORES = 0
      IF (NCO .NE. 0) NCORES = 1
      NHAM = NCORES+NPAIR+NPAIR+NSETO
C
C     ----- GENERATE KCORB -----
C
      NBASE = NCO+NOPEN
      IF(NPAIR.GT.0) THEN
         DO 340 KPAIR = 1,NPAIR
            KCORB(1,KPAIR) = NBASE+1
            KCORB(2,KPAIR) = NBASE+2
            NBASE = NBASE+2
  340    CONTINUE
      END IF
C
      IF(.NOT.OUT) RETURN
      WRITE(IW,9000) NCO,NSETO,NOPEN,NPAIR,NORB,NHAM
      WRITE(IW,9010) (NCONF(I),I=1,NORB)
      IF(NPAIR.GT.0) WRITE(IW,9020) (KCORB(1,I),I=1,NPAIR)
      IF(NPAIR.GT.0) WRITE(IW,9020) (KCORB(2,I),I=1,NPAIR)
      RETURN
C
 9000 FORMAT(1X,'NCO,NSETO,NOPEN,NPAIR,NORB=',5I5,1X,'NHAM=',I4)
 9010 FORMAT(1X,'NCONF=',15I4)
 9020 FORMAT(1X,'KCORB=',15I4)
      END
C*MODULE GVB     *DECK GVBTCH
      SUBROUTINE GVBTCH(ETOT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION H(3),E(2),V(2,2),SCR(2,8),IWRK(2)
C
      COMMON /GVBCI / CIHAM(91),CILOW(12),JPAIR,KONE,KTWO,KCORB(2,12)
C
C        ----- ADJUST TCSCF HAMILTONIAN MATRIX -----
C        THE DIAGONAL OF THE HAMILTONIAN IS MISSING THE
C        ADDITIVE CONSTANT -H00-, SO WE PUT IT BACK IN.
C        NOT ONLY IS THIS CODE GOOD JUST FOR THE TCSCF CASE, IT
C        IS NOT PARTICULARLY FIRST PRINCIPLES IN FINDING H00!
C
      H(1) = CIHAM(1)
      H(2) = CIHAM(2)
      H(3) = CIHAM(3)
      IERR=0
      CALL GLDIAG(2,2,2,H,SCR,E,V,IERR,IWRK)
      H00 = ETOT - E(1)
      CIHAM(1) = CIHAM(1) + H00
      CIHAM(3) = CIHAM(3) + H00
      RETURN
      END
C*MODULE GVB     *DECK GVBTRN
      SUBROUTINE GVBTRN(X,C,G,WRK,L0,L1,NPR,SQCDF,RSTRCT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      DIMENSION X(NPR),C(L1,L0),G(L0,L0),WRK(L0)
C
      LOGICAL RSTRCT
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
C
C     ----- APPROXIMATE ORBITAL TRANSFORMATION FOR GVB-SOSCF
C          X      = ROTATION PARAMETERS
C          G      = TRANSFORMATION MATRIX ( G = EXP(X))
C          C      = MATRIX OF MO COEFFICIENTS TO BE ROTATED
C          L0     = NUMBER OF LINEARLY INDEPENDENT BASIS FUNCS (MO'S)
C          L1     = NUMBER OF BASIS FUNCTIONS (AO'S)
C          NPR    = NUMBER OF ROTATION PARAMETERS
C
      CALL VCLR(G,1,L0*L0)
C
C     ----- FORM MATRIX G = I + X -----
C
      DO 100 I = NORB+1,L0
         NCONF(I) = NCONF(NORB)+1
  100 CONTINUE
C
      K=0
      DO 290 IHAM = 1,NHAM
         DO 280 J = 1,NORB
            IF (NCONF(J) .NE. IHAM) GO TO 280
            DO 270 I = NCO+1,L0
               IF (NCONF(I) .LE. IHAM) GO TO 270
               K=K+1
               G(J,I)= X(K)
               G(I,J)=-X(K)
  270       CONTINUE
  280    CONTINUE
  290 CONTINUE
C
      DO 310 I=1,L0
         G(I,I) = ONE
  310 CONTINUE
C
C     ----- ORTHONORMALIZE THE TRANSFORMATION MATRIX -----
C
      DO 460 I=1,L0
         DUM = DDOT(L0,G(1,I),1,G(1,I),1)
         CALL DSCAL(L0,ONE/SQRT(DUM),G(1,I),1)
         IF(I.EQ.L0) GO TO 460
         IP1=I+1
         DO 450 J=IP1,L0
            DUM = DDOT(L0,G(1,I),1,G(1,J),1)
            CALL DAXPY(L0,-DUM,G(1,I),1,G(1,J),1)
  450    CONTINUE
  460 CONTINUE
C
C     ----- POSSIBLE RESTRICTION OF ORBITAL INTERCHANGE -----
C
      IF(RSTRCT) THEN
         CALL VCLR(WRK,1,L0)
         DO 500 IHAM = 1,NHAM
            CALL SELGV1(G,WRK,L0,IHAM,RSTRCT)
  500    CONTINUE
      END IF
C
C     ----- COMPUTE MEASURE OF ORBITAL ROTATION CHANGE -----
C
      SQCDF = ZERO
      DO 640 J = 1,NORB
         JHAM = NCONF(J)
         DO 600 I = 1,L0
            SQCDF = SQCDF+G(I,J)*G(I,J)
  600    CONTINUE
         DO 620 I = 1,NORB
            IF (NCONF(I) .EQ. JHAM) SQCDF = SQCDF-G(I,J)*G(I,J)
  620    CONTINUE
  640 CONTINUE
C
C     ----- ROTATE THE ORBITALS -C- BY TRANSFORMATION -G- -----
C
      DO 750 I=1,L1
         DO 730 J=1,L0
            DUM=ZERO
            DO 720 K=1,L0
               DUM=DUM+C(I,K)*G(K,J)
  720       CONTINUE
            WRK(J)=DUM
  730    CONTINUE
         DO 740 J=1,L0
            C(I,J)=WRK(J)
  740    CONTINUE
  750 CONTINUE
      RETURN
      END
C*MODULE GVB     *DECK MPCJK
      SUBROUTINE MPCJK(TRANS,HFM,DFM,W,PTOT2,NCLF,NOPF,IFMF,L1,L2,OUT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      LOGICAL OUT
C
      DIMENSION PTOT2(NAT,16),W(N2EL),HFM(*),DFM(*),TRANS(L1,L1)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /N2ELCT/ N2EL
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      PARAMETER (HALF=0.5D+00, TWO = 2.0D+00, ONE=1.0D+00)
C
C     -----  MPCJK IS THE SEMIEMPIRICAL ANALOGUE OF HSTARG -----
C     READ MOPAC 2-EL. INTEGRALS
C
      IF(N2EL.GT.0) CALL DAREAD(IDAF,IODA,W,N2EL,53,0)
C
C     ----- CLOSED SHELL PART -----
C
      IOFF = 1
      IF(NCLF.GT.0) THEN
         CALL DENSHL(TRANS,DFM(IOFF),1,L1,L2)
C        IT IS NECESSARY TO MULTIPLY DENSITY BY 2 TO MAKE IT
C        THE SAME AS IN RHF
         CALL DSCAL(L2,TWO,DFM(IOFF),1)
         HFCO=-HALF
         CALL MPCG(HFM(IOFF),HFM(IOFF),DFM(IOFF),DUMMY,DUMMY,PTOT2,
     *             W,L2,HFCO)
         IOFF = IOFF + L2
      END IF
C
C     ----- OPEN SHELL PART -----
C
      IF(NOPF.GT.0) THEN
         ILOW = IFMF
         IHI = ILOW + NOPF - 1
         IF(NCO.EQ.0) ILOW=ILOW-1
         IF(NCO.EQ.0) IHI =IHI -1
         HFCO=ONE
         IOFF1 = IOFF
         DO 200 I = ILOW,IHI
            CALL DENSHL(TRANS,DFM(IOFF),I,L1,L2)
            CALL MPCG(HFM(IOFF1),HFM(IOFF1+L2),DFM(IOFF),DUMMY,
     *                 DUMMY,PTOT2,W,L2,HFCO)
            IOFF = IOFF + L2
            IOFF1 = IOFF1 + L2 + L2
  200    CONTINUE
      END IF
      IF(OUT) THEN
         DO 210 I=1,NCLF+NOPF
            WRITE(IW,*) 'CURRENT DENSITY MATRIX #',I
            CALL PRTRIL(DFM(1+L2*(I-1)),L1)
210      CONTINUE
      END IF
      RETURN
      END
C
C*MODULE GVB     *DECK HSTARG
      SUBROUTINE HSTARG(TRANS,HFM,DFM,XP,XK,IX,DSH,XINTS,GHONDO,
     *                  DDIJ,NIJS,NKLS,IA,NCLF,NOPF,IFMF,L1,L2,NSH2,
     *                  NINTMX,NOPK,INTTYP,MAXG,DIRSCF,SCHWRZ,
     *                  NINT,NSCHWZ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL SCHWRZ,DIRSCF,SHLOUT
C
      DIMENSION XP(NINTMX),XK(NINTMX),IX(*),NIJS(NINTMX),
     *          NKLS(NINTMX),HFM(*),DFM(*),TRANS(L1,L1),IA(L1),
     *          DSH(NSH2),XINTS(NSH2),GHONDO(MAXG),DDIJ(*)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /PCKLAB/ LABSIZ
      COMMON /SHLT  / SHLTOL,CUTOFF,ICOUNT,SHLOUT
C
      PARAMETER (HALF=0.5D+00)
C
      DATA GVB/8HGVB     /
C
C     ----- HSTARG FORMS THE GVB COULOMB AND EXCHANGE MATRICES -----
C     MULTIPLE MATRICES ARE FORMED IN ONE INTEGRAL PASS.  EACH
C     OPEN SHELL HAS ONE COULOMB, AND ONE EXCHANGE MATRIX.
C
C     IN THIS CALL TO HSTARG, NCLF (0 OR 1) CLOSED SHELL MATRICES
C     WILL BE FORMED, ALONG WITH 2*NOPF OPEN SHELL MATRICES.
C     IN THIS ROUTINE, "OPEN" MEANS OPEN SHELL OR PAIR ORBITALS.
C     THERE WILL BE NCLF+NOPF DENSITY MATRICES FORMED.
C
      IF(NCLF.EQ.0 .AND. NOPF.EQ.0) CALL ABRT
C
C     ----- ZERO OUT THE MATRICES -----
C
      NFCK = NCLF + 2*NOPF
      NDEN = NCLF + NOPF
      NZERO = NFCK*L2
      CALL VCLR(HFM,1,NZERO)
C
C     ----- LOAD UP THE DENSITY MATRICES -----
C
C     ----- CLOSED SHELL PART -----
C
      IOFF = 1
      IF(NCLF.GT.0) THEN
         CALL DENSHL(TRANS,DFM(IOFF),1,L1,L2)
         IOFF = IOFF + L2
      END IF
C
C     ----- OPEN SHELL PART -----
C
      IF(NOPF.GT.0) THEN
         ILOW = IFMF
         IHI = ILOW + NOPF - 1
         IF(NCO.EQ.0) ILOW=ILOW-1
         IF(NCO.EQ.0) IHI =IHI -1
         DO 200 I = ILOW,IHI
            CALL DENSHL(TRANS,DFM(IOFF),I,L1,L2)
            IOFF = IOFF + L2
  200    CONTINUE
      END IF
C
C     ----- DETERMINE INTEGRAL FILE TYPE -----
C
      IF(.NOT.DIRSCF) THEN
         IF(NOPK.EQ.1) GO TO 400
         IF(NOPK.NE.1) GO TO 700
      END IF
C
C     ----- DIRECT EVALUATION OF THE OPERATORS -----
C
      IF(SCHWRZ) THEN
         DUMMY = 0.0D+00
         CALL SHLDEN(GVB,DUMMY,DUMMY,DFM,DSH,IA,L1,L2,NSH2,NDEN)
      END IF
      CUTSV  = CUTOFF
      CUTOFF = MIN(CUTOFF,1.0D-10)
      CALL TWOEI(GVB,DIRSCF,.FALSE.,.FALSE.,.FALSE.,
     *           INTTYP,SCHWRZ,NINT,NSCHWZ,L1,L2,
     *           DUMMY,DUMMY,IDUMMY,1,
     *           XINTS,NSH2,GHONDO,MAXG,DDIJ,
     *           IA,DFM,HFM,DUMMY,DUMMY,DSH,DUMMY,DUMMY,1)
      CUTOFF = CUTSV
      GO TO 500
C
C     ----- PROCESS INTEGRALS NOT IN SUPERMATRIX FORM -----
C
  400 CONTINUE
      CALL SEQREW(IJK)
      NXX = 0
  410 CONTINUE
      CALL PREAD(IJK,XP,IX,NXX,NINTMX)
      IF(NXX.EQ.0) GO TO 490
      NINT = ABS(NXX)
      IF (NINT .GT. NINTMX) CALL ABRT
C
      DO 460 M = 1,NINT
C
                       NPACK = M
                       IF (LABSIZ .EQ. 2) THEN
*I32                     LABEL1 = IX( 2*NPACK - 1 )
*I32                     LABEL2 = IX( 2*NPACK     )
*I32                     IPACK = ISHFT( LABEL1, -16 )
*I32                     JPACK = IAND( LABEL1, 65535 )
*I32                     KPACK = ISHFT( LABEL2, -16 )
*I32                     LPACK = IAND( LABEL2, 65535 )
*I64                     LABEL = IX(NPACK)
*I64                     IPACK = ISHFT( LABEL, -48 )
*I64                     JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                     KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64                     LPACK = IAND( LABEL, 65535 )
                       ELSE IF (LABSIZ .EQ. 1) THEN
*I32                     LABEL = IX(NPACK)
*I32                     IPACK = ISHFT( LABEL, -24 )
*I32                     JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32                     KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I32                     LPACK = IAND( LABEL, 255 )
*I64                     IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                       LABEL = IX( NPACK/2 )
*I64                       IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64                       JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                       LPACK = IAND( LABEL, 255 )
*I64                     ELSE
*I64                       LABEL = IX( (NPACK/2)+1 )
*I64                       IPACK = ISHFT( LABEL, -56 )
*I64                       JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64                       LPACK = IAND( ISHFT( LABEL, -32 ), 255 )
*I64                     END IF
                       END IF
                       I = IPACK
                       J = JPACK
                       K = KPACK
                       L = LPACK
C
C
      NIJ = IA(I)+J
      NIK = IA(I)+K
      NIL = IA(I)+L
      NKL = IA(K)+L
      IF (J .GE. K) THEN
         NJK = IA(J)+K
         NJL = IA(J)+L
      ELSE
         NJK = IA(K)+J
         NJL = IA(J)+L
         IF (J.LT.L) NJL = IA(L)+J
      END IF
C
      VAL = XP(M)
      VAL4 = (VAL+VAL)+(VAL+VAL)
C
      IOFF1 = 0
      IOFF2 = 0
C
C     ----- CLOSED SHELL PART -----
C
      IF(NCLF.GT.0) THEN
         HFM(NIJ) = HFM(NIJ)+VAL4*DFM(NKL)
         HFM(NKL) = HFM(NKL)+VAL4*DFM(NIJ)
         HFM(NIK) = HFM(NIK)-VAL *DFM(NJL)
         HFM(NIL) = HFM(NIL)-VAL *DFM(NJK)
         HFM(NJK) = HFM(NJK)-VAL *DFM(NIL)
         HFM(NJL) = HFM(NJL)-VAL *DFM(NIK)
         IOFF1 = IOFF1+L2
         IOFF2 = IOFF2+L2
      END IF
C
C     ----- OPEN SHELL PART -----
C         THIS LOOP HAS DATA DEPENDENCIES; UNSAFE TO VECTORIZE
C
      IF(NOPF.GT.0) THEN
      DO 450 IFO = 1,NOPF
         HFM(NIJ+IOFF1)    = HFM(NIJ+IOFF1)   +VAL4*DFM(NKL+IOFF2)
         HFM(NKL+IOFF1)    = HFM(NKL+IOFF1)   +VAL4*DFM(NIJ+IOFF2)
         HFM(NIK+IOFF1+L2) = HFM(NIK+IOFF1+L2)+VAL *DFM(NJL+IOFF2)
         HFM(NIL+IOFF1+L2) = HFM(NIL+IOFF1+L2)+VAL *DFM(NJK+IOFF2)
         HFM(NJK+IOFF1+L2) = HFM(NJK+IOFF1+L2)+VAL *DFM(NIL+IOFF2)
         HFM(NJL+IOFF1+L2) = HFM(NJL+IOFF1+L2)+VAL *DFM(NIK+IOFF2)
         IOFF1 = IOFF1+L2+L2
         IOFF2 = IOFF2+L2
  450 CONTINUE
      END IF
C
  460 CONTINUE
      IF(NXX.GT.0) GO TO 410
  490 CONTINUE
      CALL SEQREW(IJK)
C
C     ----- DONE WITH INTEGRALS, DOUBLE/HALVE DIAGONALS -----
C
  500 CONTINUE
      IOFF1 = 0
      IF(NCLF.GT.0) THEN
         DO 520 M = 1,L1
            NIJ = IA(M)+M
            HFM(NIJ) = HFM(NIJ)+HFM(NIJ)
  520    CONTINUE
         IOFF1 = IOFF1+L2
      END IF
C
      IF(NOPF.GT.0) THEN
         DO 580 IFO = 1,NOPF
            DO 550 M = 2,L1
               MAX = M-1
               DO 540 N = 1,MAX
                  NIJ = IA(M)+N
                  HFM(NIJ+IOFF1) = HFM(NIJ+IOFF1)*HALF
  540          CONTINUE
  550       CONTINUE
            DO 560 M = 1,L1
               NIJ = IA(M)+M
               HFM(NIJ+IOFF1+L2) = HFM(NIJ+IOFF1+L2)+HFM(NIJ+IOFF1+L2)
  560       CONTINUE
            IOFF1 = IOFF1+L2+L2
  580    CONTINUE
      END IF
C
      RETURN
C     ******
C
C     ----- PROCESS INTEGRALS IN SUPERMATRIX FORM -----
C
C     ----- HALVE THE DIAGONAL OF THE DENSITY MATRICES -----
C
  700 CONTINUE
      IOFF = 0
      IF(NCLF.GT.0) THEN
         DO 720 J = 1,L1
            NIJ = IA(J) + J
            DFM(NIJ) = DFM(NIJ)*HALF
  720    CONTINUE
         IOFF = IOFF + L2
      END IF
C
      IF(NOPF.GT.0) THEN
         DO 780 IFO = 1,NOPF
            DO 760 J = 1,L1
               NIJ = IA(J) + J
               DFM(NIJ+IOFF) = DFM(NIJ+IOFF)*HALF
  760       CONTINUE
            IOFF = IOFF + L2
  780    CONTINUE
      END IF
C
C     ----- UNPACK A BUFFER OF LABELS -----
C
      CALL SEQREW(IJK)
      NXX = 0
*I64  MASK32 = 2**32 - 1
C
  880 CONTINUE
      CALL PKREAD(IJK,XP,XK,IX,NXX,NINTMX)
      IF (NXX .EQ. 0) GO TO 1080
      NINT = ABS(NXX)
      IF (NINT .GT. NINTMX) CALL ABRT
      DO 900 M = 1,NINT
C
                     NPACK = M
                     IF (LABSIZ .EQ. 2) THEN
*I32                   IPACK = IX( 2*NPACK - 1 )
*I32                   JPACK = IX( 2*NPACK     )
*I64                   LABEL = IX(NPACK)
*I64                   IPACK = ISHFT( LABEL, -32 )
*I64                   JPACK = IAND( LABEL, MASK32 )
                     ELSE IF (LABSIZ .EQ. 1) THEN
*I32                   LABEL = IX(NPACK)
*I32                   IPACK = ISHFT( LABEL, -16 )
*I32                   JPACK = IAND( LABEL, 65535 )
*I64                   IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                     LABEL = IX( NPACK/2 )
*I64                     IPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64                     JPACK = IAND( LABEL, 65535 )
*I64                   ELSE
*I64                     LABEL = IX( (NPACK/2)+1 )
*I64                     IPACK = ISHFT( LABEL, -48 )
*I64                     JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                   END IF
                     END IF
                     NIJS(M) = IPACK
                     NKLS(M) = JPACK
C
  900 CONTINUE
C
C     ----- BUILD THE COULOMB AND EXCHANGE MATRICES -----
C
      IOFF1 = 0
      IOFF2 = 0
C
C     ----- CLOSED SHELL PART -----
C
      IF(NCLF.GT.0) THEN
         DO 910 M = 1,NINT
            HFM(NIJS(M)) = HFM(NIJS(M)) + XP(M)*DFM(NKLS(M))
            HFM(NKLS(M)) = HFM(NKLS(M)) + XP(M)*DFM(NIJS(M))
  910    CONTINUE
         IOFF1 = IOFF1 + L2
         IOFF2 = IOFF2 + L2
      END IF
C
C     ----- OPEN SHELL PART -----
C     CONVERT FROM P AND K TO J=P+K AND K INTEGRALS
C
      IF(NOPF.GT.0) THEN
         DO 960 M = 1,NINT
            XP(M) = XP(M) + XK(M)
  960    CONTINUE
         DO 980 IFO = 1,NOPF
            CALL JKMATS(XP,XK,DFM(IOFF2+1),HFM(IOFF1+1),
     *                  HFM(IOFF1+L2+1),NIJS,NKLS,NINT)
            IOFF1 = IOFF1 + L2 + L2
            IOFF2 = IOFF2 + L2
  980    CONTINUE
      END IF
C
      IF (NXX .GT. 0) GO TO 880
C
C     ----- DONE READING PK FILE -----
C
 1080 CONTINUE
      IOFF1 = 0
      IF(NCLF.GT.0) IOFF1 = IOFF1 + L2
      IF(NOPF.GT.0) THEN
         DO 1140 IFO = 1,NOPF
            DO 1120 I = 1,L2
               HFM(I+IOFF1) = HFM(I+IOFF1)*HALF
 1120       CONTINUE
            IOFF1 = IOFF1 + L2 + L2
 1140    CONTINUE
      END IF
C
      CALL SEQREW(IJK)
      RETURN
      END
C*MODULE GVB     *DECK JKFORM
      SUBROUTINE JKFORM(JKUNIT,TRANS,HFM,DFM,XP,XK,IXPK,NIJS,NKLS,IA,
     *                  SCR,DSH,XINTS,GHONDO,DDIJ,L0,L1,L2,MXFOCK,
     *                  NPASS,NINTMX,NOPK,INTTYP,NSH2,MAXG,SCHWRZ,
     *                  DIRSCF,NINT,NSCHWZ,DBG)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL SCHWRZ,DIRSCF,DBG,OUT,GOPARR,DSKWRK,MASWRK,TDSKWRK
C
      DIMENSION TRANS(L1,L1),HFM(*),DFM(*),XP(NINTMX),XK(NINTMX),
     *          IXPK(NINTMX),NIJS(NINTMX),NKLS(NINTMX),IA(L1),
     *          SCR(L1),DSH(NSH2),XINTS(NSH2),GHONDO(MAXG),DDIJ(*)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA DBUGME/8HJKFORM  /
      DATA NONE/4HNONE/
C
C     ----- TWO ELECTRON CONTRIBUTIONS TO J AND K MATRICES -----
C     MXFOCK MATRICES CAN BE HELD IN STORAGE AT ONE TIME.  IF
C     THIS IS LESS THAN NHAM, MULTIPLE CALLS TO HSTARG, AND THUS
C     MULTIPLE READS OF THE INTEGRAL FILE, WILL BE REQUIRED.
C
      OUT = DBG  .OR.  (EXETYP.EQ.DBUGME  .AND.  MASWRK)
      TDSKWRK = DSKWRK
      IF (.NOT.DIRSCF) DSKWRK = .TRUE.
      CALL SEQREW(JKUNIT)
C
      NPASS = 1
      IFMF = 2
      NLEFT = NHAM
C
C     ----- LOOP OVER HSTARG UNTIL ALL JK MATRICES ARE BUILT -----
C
  100 CONTINUE
      NCLF = 0
      NOPF = MXFOCK
      IF(NPASS.EQ.1  .AND.  NCO.GT.0) THEN
         NCLF = 1
         NOPF = MXFOCK-1
         NLEFT = NLEFT - 1
      END IF
C
      IF(NOPF .GT. NLEFT) NOPF = NLEFT
      NLEFT = NLEFT - NOPF
C
      IF(OUT) WRITE(IW,9000) NPASS,MXFOCK,NCLF,NOPF,IFMF
 9000 FORMAT(1X,'JKFORM: NPASS,MXFOCK=',I5,I10/
     *       1X,'JKFORM: NCLF,NOPF,IFMF=',3I5)
C
C     ----- FORM THE SKELETON COULOMB AND EXCHANGE MATRICES -----
C
      IF(MPCTYP.EQ.NONE) THEN
         CALL HSTARG(TRANS,HFM,DFM,XP,XK,IXPK,DSH,XINTS,GHONDO,DDIJ,
     *               NIJS,NKLS,IA,NCLF,NOPF,IFMF,L1,L2,NSH2,
     *               NINTMX,NOPK,INTTYP,MAXG,DIRSCF,SCHWRZ,NINT,NSCHWZ)
      ELSE
         CALL MPCJK(TRANS,HFM,DFM,XP,XK,NCLF,NOPF,IFMF,L1,L2,OUT)
      ENDIF
C
      IF(OUT) THEN
         IOFF=1
         IF(NCO.GT.0) THEN
            WRITE(IW,*) 'SKELETON CORE FOCK MATRIX IN AO BASIS'
            CALL PRTRI(HFM,L1)
            IOFF=IOFF+L2
         END IF
         DO 120 I=1,NOPF
            WRITE(IW,*) 'SKELETON J MATRIX IN AO BASIS FOR OPEN SHELL',I
            CALL PRTRI(HFM(IOFF),L1)
            IOFF=IOFF+L2
            WRITE(IW,*) 'SKELETON K MATRIX IN AO BASIS FOR OPEN SHELL',I
            CALL PRTRI(HFM(IOFF),L1)
            IOFF=IOFF+L2
  120    CONTINUE
      END IF
C
C     ----- SYMMETRIZE THE J,K MATRICES -----
C
      IF (GOPARR) THEN
         CALL DDI_GSUMI(1001,NINT  ,1)
         CALL DDI_GSUMI(1002,NSCHWZ,1)
      END IF
      NFO = NCLF+2*NOPF
      IOFF = 1
      DO 230 I=1,NFO
         IF (GOPARR) CALL DDI_GSUMF(1000,HFM(IOFF),L2)
         CALL SYMH(HFM(IOFF),DFM,IA)
         IOFF = IOFF+L2
  230 CONTINUE
      DSKWRK = TDSKWRK
C
      IF(OUT) THEN
         IOFF=1
         IF(NCO.GT.0) THEN
            WRITE(IW,*) 'SYMMETRIZED CORE FOCK MATRIX IN AO BASIS'
            CALL PRTRI(HFM,L1)
            IOFF=IOFF+L2
         END IF
         DO 240 I=1,NOPF
            WRITE(IW,*) 'SYMMETRIZED 2E- J (AO BASIS) FOR OPEN SHELL',I
            CALL PRTRI(HFM(IOFF),L1)
            IOFF=IOFF+L2
            WRITE(IW,*) 'SYMMETRIZED 2E- K (AO BASIS) FOR OPEN SHELL',I
            CALL PRTRI(HFM(IOFF),L1)
            IOFF=IOFF+L2
  240    CONTINUE
      END IF
C
C     ----- TRANSFORM TO THE MO BASIS -----
C
      IOFF = 0
      IF (NPASS .GT. 1) GO TO 360
C
C     ----- CLOSED SHELL PART -----
C     CORE HAMILTONIAN INTEGRALS ARE READ, AND ADDED TO THE
C     TWO ELECTRON CONTRIBUTION TO THE J,K MATRICES.
C
C     IF THERE IS NO CORE, THE FIRST PART OF HFM IS AN OPEN
C     SHELL J, AND MUST BE SAVED FOR LATER.  IN THIS CASE,
C     ONLY THE ONE ELECTRON INTEGRALS CONTRIBUTE TO THE CORE.
C
      IF(NCO.GT.0) THEN
         IOFF = IOFF + L2
         CALL DAREAD(IDAF,IODA,DFM,L2,11,0)
         CALL VADD(HFM,1,DFM,1,HFM,1,L2)
      ELSE
         CALL DCOPY(L2,HFM,1,DFM(L2+1),1)
         CALL DAREAD(IDAF,IODA,HFM,L2,11,0)
      END IF
C
      IF(OUT) THEN
         WRITE(IW,*) 'COMPLETE CORE FOCK MATRIX IN AO BASIS'
         CALL PRTRI(HFM,L1)
      END IF
      IF (.NOT.DIRSCF) DSKWRK = .TRUE.
      CALL TFTRI(DFM,HFM,TRANS,SCR,L0,L1,L1)
      IF(OUT) THEN
         WRITE(IW,*) 'COMPLETE CORE FOCK MATRIX IN MO BASIS'
         CALL PRTRI(DFM,L1)
      END IF
      CALL SQWRIT(JKUNIT,DFM,L2)
      IF(NCO.EQ.0) CALL DCOPY(L2,DFM(L2+1),1,HFM,1)
C
C     ----- OPEN SHELL PART -----
C
  360 CONTINUE
      IF (.NOT.DIRSCF) DSKWRK = .TRUE.
      IF(NOPF.GT.0) THEN
         DO 380 IFO = 1,NOPF
            CALL TFTRI(DFM,HFM(1+IOFF),TRANS,SCR,L0,L1,L1)
            IF(OUT) THEN
               WRITE(IW,*) 'J OPERATOR IN MO BASIS FOR SHELL',IFO
               CALL PRTRI(DFM,L1)
            END IF
            CALL SQWRIT(JKUNIT,DFM,L2)
            CALL TFTRI(DFM,HFM(1+IOFF+L2),TRANS,SCR,L0,L1,L1)
            IF(OUT) THEN
               WRITE(IW,*) 'K OPERATOR IN MO BASIS FOR SHELL',IFO
               CALL PRTRI(DFM,L1)
            END IF
            CALL SQWRIT(JKUNIT,DFM,L2)
            IOFF = IOFF + L2 + L2
  380    CONTINUE
      END IF
      DSKWRK = TDSKWRK
C
C     ----- MAKE ANOTHER INTEGRAL PASS IF NECESSARY -----
C
      IF(NLEFT.EQ.0) RETURN
      NPASS = NPASS + 1
      IFMF = IFMF + NOPF
      GO TO 100
      END
C*MODULE GVB     *DECK JKMATS
      SUBROUTINE JKMATS(XP,XK,D,H1,H2,NIJS,NKLS,NINT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION XP(NINT),XK(NINT),NIJS(NINT),NKLS(NINT)
      DIMENSION D(*),H1(*),H2(*)
C
C     -----  FORM COULOMB AND EXCHANGE OPERATORS -----
C     H1 IS J, H2 IS K.  J AND K INTEGRALS ARE IN XP AND XK.
C
      DO 100 M = 1,NINT
         NIJ = NIJS(M)
         NKL = NKLS(M)
         H1(NIJ) = H1(NIJ)+XP(M)*D(NKL)
         H2(NIJ) = H2(NIJ)+XK(M)*D(NKL)
         H1(NKL) = H1(NKL)+XP(M)*D(NIJ)
         H2(NKL) = H2(NKL)+XK(M)*D(NIJ)
  100 CONTINUE
      RETURN
      END
C*MODULE GVB     *DECK LVLGVB
      SUBROUTINE LVLGVB(H,NOCC,NHOCC,N,DE,DEP,ITERL,NCALL)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION H(*)
C
      LOGICAL VSHIFT,EXTRAH
C
      SAVE DIFFP,PSHIFT
C
      COMMON /ACONV / RRSHFT,EXTTOL,DMPTOL,VSHTOL,IEXTIN
      COMMON /CONV  / ACURCY,EN,ETOT,EHF,EHF0,DIFF,ITER,ICALP,ICBET
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
      PARAMETER (ZERO=0.0D+00,  HOCCF=0.5D+00)
      PARAMETER (RBASE=2.0D+00, RFACT=0.4D+00, NSAME=4)
      PARAMETER (EFACT=5.0D+00, RBADE=0.8D+00, BADEF=2.0D+00)
      PARAMETER (DRAT=2.0D+00,  RBADD=0.2D+00, BADDF=2.0D+00)
      PARAMETER (RDROP=4.0D+00, RMIN=0.1D+00,  RMAX=10.0D+00)
C
      DATA DIFFP,PSHIFT/0.0D+00,0.0D+00/
C
C     ----- SHIFT THE DIAGONAL ELEMENTS OF THE FOCK MATRIX IN
C           THE MOLECULAR BASIS SET BY THE CHANGE IN
C           THE ENERGY AND DENSITY MATRIX.
C
      IF (ITER .EQ. 1) GO TO 240
      IF (NCALL .NE. 1) GO TO 140
      ITERL = ITERL+1
      EXTRAH = MOD(MCONV,2) .EQ. 1
      VSHIFT = MOD(MCONV,8) .GE. 4
C
C     ----- DECREASE THE LEVEL AFTER NSAME ITERATIONS -----
C
      PWR = ITER-2
      PWR = PWR/NSAME
      RSHC = RBASE*RFACT**PWR
C
C     ----- USE ENERGY AS A CRITERION TO DETERMINE THE LEVEL -----
C
      RSHE = MIN(RBASE,ABS(DE)*EFACT)
      IF (DE .GT. ZERO) RSHE = MAX(RSHE,RBADE,RRSHFT*BADEF)
      IF (DEP .GT. ZERO .AND. DE+DEP .GT. ZERO) RSHE = MAX(RSHE,
     +     RRSHFT*BADEF)
C
C     ----- USE THE DENSITY MATRIX DIFFERENCE AS A CRITERION -----
C
      RSHD = ZERO
      IF (DIFF .GT. DRAT*DIFFP .AND. RRSHFT .EQ. PSHIFT) RSHD =
     +     MAX(RBADD,RRSHFT*BADDF)
C
C     ----- TAKE THE HIGHEST LEVEL SHIFT -----
C
      RHIGH = MAX(RSHC,RSHE,RSHD)
      PSHIFT = RRSHFT
C
C     ----- LET THE LEVEL DROP AFTER NSAME ITERATIONS -----
C
      IF (MOD(ITERL,NSAME) .EQ. 0) RRSHFT = RRSHFT/RDROP
      RRSHFT = MAX(RHIGH,RRSHFT)
C
C     ----- BOUND THE LEVEL BETWEEN LIMITS -----
C
      IF (RRSHFT .LE. RMIN) RRSHFT = RMIN
      IF (RRSHFT .GT. RMAX) RRSHFT = RMAX
C
C     ----- CHECK FOR CLOSE TO CONVERGENCE OR AUTOMATIC
C           SWITCHING TO EXTRAPOLATION -----
C
      IF (DIFF.LT.ACURCY  .AND.  ITERL.NE.0) RRSHFT = ZERO
      IF (.NOT.VSHIFT) RRSHFT = ZERO
      IF (PSHIFT .EQ. ZERO) GO TO 100
      IF (EXTRAH .AND. RRSHFT .LT. VSHTOL .AND. ABS(DE) .LT. DMPTOL)
     +     RRSHFT = ZERO
      GO TO 120
C
  100 IF (EXTRAH .AND. ABS(DE) .LT. EXTTOL) RRSHFT = ZERO
C
  120 IF (RRSHFT .NE. PSHIFT) ITERL = 0
C
C     ----- ADD THE LEVEL SHIFT TO THE FOCK MATRIX -----
C
  140 CONTINUE
      IF (NHOCC .EQ. 0) GO TO 180
C
C     ----- SHIFT ANY HALF FILLED ORBITALS BY RRSHFT*HOCCF -----
C
      IS = NOCC+1
      IE = NHOCC
      DO 160 I = IS,IE
      II = I*(I+1)/2
      H(II) = H(II)+RRSHFT*HOCCF
  160 CONTINUE
C
C     ----- SHIFT ANY VIRTUAL ORBITALS BY RRSHFT -----
C
  180 IS = NOCC+NHOCC+1
      IE = N
      IF (IS .GT. IE) GO TO 220
      DO 200 I = IS,IE
      II = I*(I+1)/2
      H(II) = H(II)+RRSHFT
  200 CONTINUE
  220 IF (NCALL .EQ. 1) DIFFP = DIFF
      RETURN
C
  240 ITERL = 1
      PSHIFT = RRSHFT
      DIFFP = DIFF
      RETURN
      END
C*MODULE GVB     *DECK NEWCI
      SUBROUTINE NEWCI(H,DJ,DK,L1,CONCI,HCI,EIG,ECI,SCR,IWRK,
     *                 IA,NTERM,OUT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL OUT,GOPARR,DSKWRK,MASWRK
C
      DIMENSION DJ(*),DK(*),H(*),HCI(3),EIG(NTERM,2),ECI(2),
     *          SCR(NTERM,8),IA(L1),IWRK(L1)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /GVBCI / CIHAM(91),CILOW(12),JPAIR,KONE,KTWO,KCORB(2,12)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00, THRESH=1.0D-08)
C
C     ----- CALCULATE A NEW SET OF CI COEFFICIENTS FOR THE
C           JPAIR'TH PAIR -----
C
      IPICK = 1
      ICORE = 1
      IF(NCO.EQ.0) ICORE=0
C
      I1 = IA(KONE)+KONE
      I3 = IA(KTWO)+KTWO
C
C     ----- INITIALIZE CI MATRIX WITH HCORE -----
C
      NN = (NTERM*NTERM+NTERM)/2
      CALL VCLR(HCI,1,NN)
      HCI(1) = TWO*H(I1)
      HCI(2) = ZERO
      HCI(3) = TWO*H(I3)
C
C      ----- PROCESS J AND K MATRICES -----
C
      DO 120 II = 1,NORB
         ISH = NCONF(II)
         IX = IA(KONE)+II
         IF (II .GT. KONE) IX = IA(II)+KONE
         IY = IA(KTWO)+II
         IF (II .GT. KTWO) IY = IA(II)+KTWO
         IF (ISH .EQ. ICORE) GO TO 120
         IF (II.EQ.KONE  .OR.  II.EQ.KTWO) GO TO 120
C
C        ----- ADD POTENTIAL TERMS FROM OTHER SHELLS -----
C
         CNST = F(ISH)
         HCI(1) = HCI(1)+CNST*TWO*(TWO*DJ(IX)-DK(IX))
         HCI(3) = HCI(3)+CNST*TWO*(TWO*DJ(IY)-DK(IY))
  120 CONTINUE
      IX = IA(KONE)+KONE
      IY = IA(KTWO)+KTWO
      IZ = IA(KONE)+KTWO
      IF (KTWO .GT. KONE) IZ = IA(KTWO)+KONE
      HCI(1) = HCI(1)+DJ(IX)
      HCI(2) =        DK(IZ)
      HCI(3) = HCI(3)+DJ(IY)
      HSTR = MIN(HCI(1),HCI(3))
C
      IF(OUT) THEN
         IB = 0
         WRITE (IW,9008)
         DO 140 J = 1,2
            IAA = IB+1
            IB = IB+J
            WRITE (IW,9028) (HCI(I),I = IAA,IB)
  140    CONTINUE
      END IF
C
C     ----- SAVE HAMILTONIAN ELEMENTS -----
C     ONLY 2X2 BLOCKS ARE COMPUTED HERE, SO WE STORE ONLY
C     A SUBBLOCK OF H.  THE TRUE HAMILTONIAN SHOULD BE OF
C     2**NPAIR X 2**NPAIR, BUT HERE THE DIMENSION IS JUST
C     NPAIR+1 X NPAIR+1, AND ONLY THE FIRST COLUMN AND THE
C     DIAGONAL OF THIS SUBMATRIX IS COMPUTED.   THE NON-PAIR
C     CONTRIBUTION (H00) IS NOT INCLUDED AS IT DOES NOT HAVE
C     ANY EFFECT ON THE EIGENVECTORS WE ARE LOOKING FOR.
C
      I = 1
      J = JPAIR + 1
      II = IA(I) + I
      IJ = IA(J) + I
      JJ = IA(J) + J
      CIHAM(II) = HCI(1)
      CIHAM(IJ) = HCI(2)
      CIHAM(JJ) = HCI(3)
C
C     ----- GET NEW CI COEFFICIENTS -----
C
      IERR = 0
      CALL GLDIAG(NTERM,NTERM,NTERM,HCI,SCR,ECI,EIG,IERR,IWRK)
      IF (IERR .NE. 0) CALL ABRT
      SIGNSW = EIG(1,IPICK)
      DO 180 K = 1,NTERM
         IF (SIGNSW .LT. ZERO) EIG(K,IPICK) = -EIG(K,IPICK)
         CONCI = CONCI+(CICOEF(K,JPAIR)-EIG(K,IPICK))**2
         CICOEF(K,JPAIR) = EIG(K,IPICK)
  180 CONTINUE
      CILOW(JPAIR) = ECI(1)-HSTR
      IF(OUT) THEN
         WRITE (IW,9048)
         DO 200 J = 1,NTERM
            WRITE (IW,9068) J,ECI(J),(EIG(I,J),I = 1,NTERM)
  200    CONTINUE
      END IF
      LONE = NCONF(KONE)
      LTWO = NCONF(KTWO)
      F(LONE) = EIG(1,IPICK)**2
      F(LTWO) = EIG(2,IPICK)**2
      IF (F(LONE) .LT. THRESH) GO TO 240
      IF (F(LTWO) .LT. THRESH) GO TO 240
      IF ((EIG(1,IPICK)/EIG(2,IPICK)) .GT. ZERO) GO TO 260
      RETURN
C
  240 CONTINUE
      IF (MASWRK) WRITE (IW,9088) THRESH,JPAIR,F(LONE),F(LTWO)
      CALL ABRT
      STOP
C
  260 CONTINUE
      IF (MASWRK) WRITE (IW,9108) (EIG(K,IPICK),K = 1,NTERM)
      CALL ABRT
      STOP
C
 9008 FORMAT(/5X,'CI MATRIX')
 9028 FORMAT(5X,1P,5E15.8)
 9048 FORMAT(//5X,12HCI SOLUTIONS  /)
 9068 FORMAT(5X,I5,5X,8HENERGY= ,F20.10/(5X,5E15.8))
 9088 FORMAT(//5X,'JOB ABORTED WITH OCCUPATION NUMBER BELOW ',1P,E12.3/
     *         5X,'PAIR',I4,' HAS OCC NOS=',1P,2E15.8)
 9108 FORMAT(//5X,'CI COEFS ARE SAME SIGN - GOODBYE' ,2E15.8//)
      END
C*MODULE GVB     *DECK OPTCI
      SUBROUTINE OPTCI(JKUNIT,H,DJ,SCR,IWRK,IA,L1,L2,OUT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL OUT
C
      DIMENSION DJ(*),H(*),SCR(L1,8),IA(L1),IWRK(L1)
      DIMENSION HCI(3),EIG(2,2),ECI(2)
C
      COMMON /CONV  / DENTOL,EN,ETOT,EHF,EHFO,DIFF,ITER,ICALP,ICBET
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /GVBCI / CIHAM(91),CILOW(12),JPAIR,KONE,KTWO,KCORB(2,12)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
      PARAMETER (ZERO=0.0D+00, PT01=1.0D-02)
C
C     ----- LOOP OVER ALL PAIRS AND OPTIMIZE THE CI COEFFICIENTS -----
C
      IF (NPAIR .EQ. 0) RETURN
      ICORE=1
      IF(NCO.EQ.0) ICORE=0
      TEST = DENTOL*PT01
      NTERM = 2
      IAA = L2+1
C
C     ----- GET J(I,J) AND K(I,J) INTEGRALS IN PLACE -----
C
      CALL SEQREW(JKUNIT)
      CALL SQREAD(JKUNIT,H,L2)
      IOLD = 0
      DO 140 I = 1,NORB
         IF (NCONF(I).EQ.ICORE) GO TO 140
         IF (NCONF(I).NE.IOLD) THEN
            CALL SQREAD(JKUNIT,H,L2)
            CALL SQREAD(JKUNIT,H(IAA),L2)
         END IF
         DO 120 J = 1,NORB
            IF (NCONF(J).EQ.ICORE) GO TO 120
            JJ = IA(J)+J
            KK = JJ+L2
            IJ = IA(I)+J
            IF (J .GT. I) IJ = IA(J)+I
            DJ(IJ   ) = H(JJ)
            DJ(IJ+L2) = H(KK)
  120    CONTINUE
         IOLD = NCONF(I)
  140 CONTINUE
      CALL SEQREW(JKUNIT)
      CALL SQREAD(JKUNIT,H,L2)
C
C     ----- START CI ITERATION OVER PAIRS -----
C
      ITERCI = 0
  160 CONTINUE
      ITERCI = ITERCI+1
      CONCI = ZERO
C
C       ----- OPTIMIZE THE JPAIR'TH PAIR -----
C
      DO 180 JPAIR = 1,NPAIR
         KONE = KCORB(1,JPAIR)
         KTWO = KCORB(2,JPAIR)
         CALL NEWCI(H,DJ,DJ(IAA),L1,CONCI,HCI,EIG,ECI,
     *              SCR,IWRK,IA,NTERM,OUT)
  180 CONTINUE
      IF (NPAIR .EQ. 1) CONCI = ZERO
C
C     ----- TEST FOR CONVERGENCE OF THE CI OPTIMIZATION -----
C
      IF (CONCI .GT. TEST) GO TO 160
C
C     ----- NORMALIZE THE CI COEFFICIENTS AND UPDATE THE F, ALPHA,
C           AND BETA MATRICES -----
C
      CALL CIEXPR(KCORB,IA,L1)
      IF(.NOT.OUT) RETURN
C
      WRITE (IW,9008) ITERCI
      WRITE (IW,9028)
      DO 200 J = 1,NHAM
      WRITE (IW,9048) J,F(J)
  200 CONTINUE
      WRITE(IW,9058)
      CALL PRTRI(ALPHA,NHAM)
      CALL PRTRI(BETA,NHAM)
      RETURN
C
 9008 FORMAT(5X,21HCI PART CONVERGED IN  ,I3,2X,10HITERATIONS  )
 9028 FORMAT(5X,12HENERGY COEFS  )
 9048 FORMAT(5X,6HSHELL ,I2,2X,3HF= ,F15.8)
 9058 FORMAT(1X,'UPDATED ALPHA AND BETA MATRICES')
      END
C*MODULE GVB     *DECK REDDEN
      SUBROUTINE REDDEN(TRANS,DA,DB,L1,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
      DIMENSION TRANS(L1,1),DA(L2),DB(L2)
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00)
C
C     ----- CALCULATE THE REDUCED DENSITY MATRIX -----
C
      CALL VCLR(DA,1,L2)
      CALL VCLR(DB,1,L2)
C
C     ----- GET ALPHA PART FIRST -----
C
      IJ = 0
      DO 280 I = 1,L1
         DO 270 J = 1,I
            IJ = IJ + 1
            SUM = ZERO
            DO 260 K = 1,NORB
               SUM = SUM+F(NCONF(K))*TRANS(I,K)*TRANS(J,K)
  260       CONTINUE
            DA(IJ) = SUM
  270    CONTINUE
  280 CONTINUE
C
C     ----- ALPHA=BETA FOR NO OPEN SHELLS, OR SINGLET STATE -----
C
      IF(NOPEN.EQ.0  .OR.  (NOPEN.EQ.2  .AND.  MUL.EQ.1)) THEN
         CALL DCOPY(L2,DA,1,DB,1)
         GO TO 400
      END IF
C
C     ----- IF THERE ARE OPEN SHELLS, FIRST GET THE OPEN PART -----
C     NOTE THAT THIS ROUTINE ASSUMES A HIGH SPIN COUPLING OF THE
C     OPEN SHELLS WHOSE FRACTIONAL OCCUPANCY IS IN THE F ARRAY.
C     THE GVB CODE WORKS WITH THE TOTAL DENSITY, SO THE ENERGY FOR
C     LOW SPIN COUPLING WILL BE CORRECT.  LATER PROPERTIES WON'T BE.
C     FOR EXAMPLE, A DOUBLET WITH 3 UNPAIRED ELECTRONS HAS WRONG SPIN
C     DENSITIES COMPUTED LATER ON.
C
      ILO = NCO+1
      IHI = NCO+NOPEN
      IJ = 0
      DO 320 I = 1,L1
         DO 310 J = 1,I
            IJ = IJ + 1
            SUM = ZERO
            DO 300 K = ILO,IHI
               SUM = SUM+F(NCONF(K))*TRANS(I,K)*TRANS(J,K)
  300       CONTINUE
            DB(IJ) = SUM
  310    CONTINUE
  320 CONTINUE
C
C     ----- ADD THIS TO ALPHA TO GET CORRECT ALPHA, AND THEN
C     SUBTRACT TWICE THE OPEN SHELL FROM ALPHA TO GET BETA -----
C
      CALL VADD(DA,1,DB,1,DA,1,L2)
      CALL DSCAL(L2,-TWO,DB,1)
      CALL VADD(DA,1,DB,1,DB,1,L2)
C
C     ----- SAVE THE REDUCED DENSITY MATRICES -----
C
  400 CONTINUE
      CALL DAWRIT(IDAF,IODA,DA,L2,16,0)
      CALL DAWRIT(IDAF,IODA,DB,L2,20,0)
      RETURN
      END
C*MODULE GVB     *DECK RHFGVB
      SUBROUTINE RHFGVB
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH,
     *           MXATM=2000, MXAO=8192)
C
      LOGICAL PACK2E,OUT,SOME,DODIIS,SCHWRZ,DIRSCF,FDIFF,
     *        GOPARR,DSKWRK,MASWRK,SOSCF
C
      DIMENSION IA(MXAO)
C
      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /N2ELCT/ N2EL
      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 /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA DEBUG/8HDEBUG   /, DBGME/8HGVB     /
      DATA NONE/4HNONE/
C
C     ----- MEMORY ALLOCATION FOR ROHF-GVB -----
C
      OUT = NPRINT.EQ.5  .OR.  EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBGME
     *      .AND. MASWRK
      SOME = NPRINT.NE.-5 .AND. MASWRK
      IF (MASWRK) WRITE (IW,9000)
C
C     ----- SET UP TRIANGULAR ADDRESS ARRAY -----
C
      DO 100 I = 1,MXAO
         IA(I) = (I*I-I)/2
  100 CONTINUE
C
C     ----- SET UP LENGTH VALUES -----
C
      L0 = NQMT
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
C
      SCHWRZ = ISCHWZ.EQ.1
      IF(DIRSCF) THEN
         CALL BASCHK(LMAX)
                       NANGM =  4
         IF(LMAX.EQ.2) NANGM =  6
         IF(LMAX.EQ.3) NANGM = 10
         IF(LMAX.EQ.4) NANGM = 15
         IF(LMAX.EQ.5) NANGM = 21
         IF(LMAX.EQ.6) NANGM = 28
         MAXG = NANGM**4
         NSH2=1
         IF(SCHWRZ) NSH2 = (NSHELL*NSHELL+NSHELL)/2
         MXINT=1
      ELSE
         MAXG=1
         NSH2=1
         MXINT=NINTMX
      END IF
C
C     ----- SET UP ROHF-GVB DATA IN /GVBWFN/  -----
C     SOME OF THIS IS ALREADY DONE IN ROUTINE SCFIN
C
      CALL GVBSET(OUT)
C
C      ----- ALLOCATE FAST MEMORY -----
C
      CALL VALFM(LOADFM)
      LTRANS = LOADFM + 1
      LSCR   = LTRANS + L3
      LIWRK  = LSCR   + L1*8
      LE     = LIWRK  + L1
      LEIG   = LE     + L1
      LS     = LEIG   + L1
      LQ     = LS     + L2
      LDSH   = LQ     + L3
      LXINTS = LDSH   + NSH2
      LGHOND = LXINTS + NSH2
      LDDIJ  = LGHOND + MAXG
      IF(MPCTYP.EQ.NONE) THEN
         LIXPK  = LDDIJ  + 49*MXG2
         LXP    = LIXPK  + MXINT
         LXK    = LXP    + MXINT
         LNIJS  = LXK    + MXINT
         LNKLS  = LNIJS  + MXINT
         LWRK1  = LNKLS  + MXINT
      ELSE
C             XP HOLDS W (2-EL. INT.), XK- BUFFER FOR PTOT2
C             LIXPK,LNIJS,LMKLS,LWRK1 ARE NOT USED
         LXP  =LDDIJ+1
         LXK  =LXP+N2EL
         LWRK1=LXK+NAT*16
         LIXPK=0
         LNIJS=0
         LNKLS=0
      END IF
      LWRK2  = LWRK1  + L3
      LWRK3  = LWRK2  + L3
      LAST   = LWRK3  + L3
C
C     GET MEMORY FOR DIIS
C
      DODIIS = MOD(MCONV,32) .GE. 16
      IF(NPAIR.LE.1  .AND.  DODIIS) THEN
         MAXIT1 = MAXIT
      ELSE
         MAXIT1 = 1
      END IF
      MAXIT2 = (MAXIT1*MAXIT1+MAXIT1)/2
      MXDII3 = MAXDII*MAXDII
      MAXIO  = 4*MAXIT1
C
      LADIIS = LAST
      LXDIIS = LADIIS + MXDII3
      LPDIIS = LXDIIS + MAXIT1
      LBDIIS = LPDIIS + MAXIT1
      LIODII = LBDIIS + MAXIT2
      LAST   = LIODII + MAXIO
C
C     ----- GET MEMORY FOR SECOND ORDER SCF -----
C
      NPR=NORB*(L0-NORB)
      DO 200 IHAM = 1,NHAM
         DO 190 J = 1,NORB
            IF (NCONF(J) .NE. IHAM) GO TO 190
            DO 180 I = NCO+1,NORB
               IF (NCONF(I) .LE. IHAM) GO TO 180
                  NPR=NPR+1
  180       CONTINUE
  190    CONTINUE
  200 CONTINUE
C
      SOSCF  = MOD(MCONV,128).GE. 64
      IF(SOSCF) THEN
         LHESS  = LAST
         LGRAD  = LHESS  + NPR
         LPGRAD = LGRAD  + NPR
         LDISPL = LPGRAD + NPR
         LWRK4  = LDISPL + NPR
         LWRK5  = LWRK4  + NPR
         LWRK6  = LWRK5  + NPR
         LDIAGH = LWRK6  + NPR
         LDIAGJ = LDIAGH + NHAM*L1
         LDIAGK = LDIAGJ + NHAM*L1
         LAST   = LDIAGK + NHAM*L1
      ELSE
         LHESS  = LAST
         LGRAD  = LAST
         LPGRAD = LAST
         LDISPL = LAST
         LWRK4  = LAST
         LWRK5  = LAST
         LWRK6  = LAST
         LDIAGH = LAST
         LDIAGJ = LAST
         LDIAGK = LAST
         NPR = 1
      END IF
C
      NEED1 = LAST - LOADFM - 1
      CALL GETFM(NEED1)
C
C     ----- ALLOCATE FOCK (COULOMB/EXCHANGE) AND DENSITY MATRICES -----
C
      CALL VALFM(LOADFM)
      CALL GOTFM(NGOT)
      MXL2 = NGOT/L2
      MXD  = MXL2/3
      MXF  = MXD+MXD
      NCORES = 0
      IF(NCO.GT.0) NCORES=1
      NFOCK  = NCORES + 2*(NHAM-NCORES)
      MXDENS = MIN( NHAM*L2,MXD*L2)
      MXFOCK = MIN(NFOCK*L2,MXF*L2)
      MXD    = MAX(MXDENS,3*L3)
C
      IF(DIRSCF  .AND.  NFOCK.GT.MXF) THEN
         NEEDX = 3*L2*(NFOCK-MXF)
         IF (MASWRK) WRITE(IW,9020) NGOT,MXF,NFOCK,NEEDX
         CALL ABRT
         STOP
      END IF
C
C        NOTE THAT DK AND H2 ARE PARTS OF DJ AND H, RESPECTIVELY.
C HUI    GIVE LH2 SOME LENGTH FOR SOLVATION MODELS
C        OR LAST=LH2 WHEN MXFOCK=L2
C
      LDJ    = LOADFM + 1
      LDK    = LDJ    + L2
      LH     = LDJ    + MXD
      LH2    = LH     + L2
      LAST   = LH2    + MXFOCK
      NEED2  = LAST   - LOADFM - 1
      CALL GETFM(NEED2)
C
      NEED = NEED1 + NEED2
      IF(SOME) WRITE(IW,9010) NEED
C
      MXFOCK = MXDENS/L2
C
      CALL GVBITR(X(LTRANS),X(LDJ),X(LDK),X(LH),X(LH2),X(LE),
     *            X(LEIG),X(LXP),X(LXK),X(LIXPK),X(LNIJS),X(LNKLS),
     *            X(LS),X(LQ),X(LSCR),X(LIWRK),IA,X(LDSH),X(LXINTS),
     *            X(LGHOND),X(LDDIJ),
     *            MXFOCK,NINTMX,L0,L1,L2,L3,NSH2,MAXG,INTTYP,
     *            DIRSCF,SCHWRZ,X(LWRK1),X(LWRK2),X(LWRK3),
     *            X(LADIIS),X(LXDIIS),X(LPDIIS),X(LBDIIS),X(LIODII),
     *            X(LHESS),X(LGRAD),X(LPGRAD),X(LDISPL),X(LWRK4),
     *            X(LWRK5),X(LWRK6),X(LDIAGH),X(LDIAGJ),X(LDIAGK),
     *            NHAM,
     *            MAXIT1,MAXIT2,MXDII3,MAXIO,NPR)
C
C     ----- RESET MEMORY -----
C
      CALL RETFM(NEED2)
      CALL RETFM(NEED1)
      IF (MASWRK)
     * WRITE(IW,FMT='(/'' ... END OF ROHF-GVB SCF CALCULATION ...'')')
      CALL TIMIT(1)
      RETURN
C
 9000 FORMAT(/10X,24(1H-)/10X,'ROHF-GVB SCF CALCULATION'/10X,24(1H-))
 9010 FORMAT(1X,'GVB STEP WILL USE',I10,' WORDS OF MEMORY.')
 9020 FORMAT(1X,'DIRECT GVB MUST HOLD ALL FOCK AND DENSITY MATRICES'/
     *       1X,'IN MEMORY.  GAMESS HAS',I10,' WORDS OF FREE MEMORY,'/
     *       1X,'WHICH WILL HOLD JUST',I5,' OF THE ',I5,' MATRICES.'/
     *       1X,'ADD AT LEAST',I10,' WORDS MORE MEMORY AND RESUBMIT.')
      END
C*MODULE GVB     *DECK ROTB
      SUBROUTINE ROTB(JKUNIT,HOLD,H,H2,ROTWT,IA,L1,L2,OUT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL OUT
C
      DIMENSION HOLD(L1,L1),H(L2),H2(L2),IA(L1)
C
      COMMON /CONV  / DENTOL,EN,ETOT,EHF,EHFO,DIFF,ITER,ICALP,ICBET
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
      PARAMETER (ZERO=0.0D+00, HALF=0.5D+00, TWO=2.00D+00, TEN=1.0D+01)
      PARAMETER (TM9=1.0D-09)
C
C     ----- ROTATE THE OCCUPIED ORBITALS AMONG THEMSELVES -----
C
      SQCDF = ZERO
      ROTTOL = DENTOL
      IF (ROTTOL .GT. TM9) ROTTOL = TM9
      ATEST = ROTTOL/TEN
      DTEST = ATEST/TEN
      ROTMAX = HALF
      CALL VCLR(HOLD,1,L1*L1)
C
      CALL SEQREW(JKUNIT)
      CALL SEQADV(JKUNIT)
C
C     ----- CALCULATE B(IJ,0) -----
C
      KCOLD = 0
      DO 220 I = 1,NORB
         KC = NCONF(I)
         IF (NCO.GT.0  .AND.  KC.EQ.1) GO TO 220
         IF (KC.NE.KCOLD) CALL SQREAD(JKUNIT,H,L2)
         DO 160 J = 1,I
            LC = NCONF(J)
            KCLC = IA(KC)+LC
            KCKC = IA(KC)+KC
            LCLC = IA(LC)+LC
            IF (KC .EQ. LC) GO TO 160
            CJ = BETA(KCKC)-BETA(KCLC)+BETA(LCLC)-BETA(KCLC)
            IJ = IA(J)+J
            HOLD(J,I) = HOLD(J,I)+CJ*H(IJ)
  160    CONTINUE
         IF (KC.NE.KCOLD) CALL SQREAD(JKUNIT,H2,L2)
         DO 200 J = 1,I
            LC = NCONF(J)
            KCLC = IA(KC)+LC
            KCKC = IA(KC)+KC
            LCLC = IA(LC)+LC
            IF (KC .EQ. LC) GO TO 200
            CK = BETA(KCKC)-BETA(KCLC)+TWO*ALPHA(KCKC)-TWO*ALPHA(KCLC)+
     +           BETA(LCLC)-BETA(KCLC)+TWO*ALPHA(LCLC)-TWO*ALPHA(KCLC)
            IJ = IA(J)+J
            HOLD(J,I) = HOLD(J,I)+CK*H2(IJ)
  200    CONTINUE
         KCOLD = KC
  220 CONTINUE
C
C     ----- CALCULATE B(IJ,IJ) -----
C
      DO 340 IHAM = 1,NHAM
         CALL SQREAD(JKUNIT,H,L2)
         DO 300 I = 1,NORB
            NX = NCONF(I)
            IF (NX .NE. IHAM) GO TO 300
            DO 280 J = 1,NORB
               I1 = IA(I)+I
               I2 = IA(J)+J
               IF(I.GE.J) THEN
                  IJ = IA(I)+J
                  HOLD(I,J) = HOLD(I,J)-F(NX)*H(IJ)
                  HOLD(J,I) = HOLD(J,I)+F(NX)*(H(I2)-H(I1))*ROTWT
               ELSE
                  IJ = IA(J)+I
                  HOLD(J,I) = HOLD(J,I)+F(NX)*H(IJ)
                  HOLD(I,J) = HOLD(I,J)+F(NX)*(H(I2)-H(I1))*ROTWT
               END IF
  280       CONTINUE
            HOLD(I,I) = ZERO
  300    CONTINUE
  340 CONTINUE
C
C     ----- PERFORM THE ACTUAL ROTATION -----
C
      DO 400 J = 1,NORB
         JC = NCONF(J)
         DO 380 I = 1,J
            IC = NCONF(I)
            IF (JC .EQ. IC) GO TO 380
            ANUM = HOLD(J,I)
            DENOM = ABS(HOLD(I,J))
            IF (ABS(DENOM) .LT. DTEST) GO TO 360
            IF (ABS(ANUM) .LT. ATEST) GO TO 360
            THETA = -ANUM/DENOM
            IF ( ABS(THETA) .GT. ROTMAX)
     *           THETA = THETA/ ABS(THETA)*ROTMAX
            HOLD(J,I) = THETA
            HOLD(I,J) = -HOLD(J,I)
            SQCDF = SQCDF+4.00D+00*HOLD(J,I)*HOLD(J,I)
            GO TO 380
C
  360       CONTINUE
            HOLD(J,I) = ZERO
            HOLD(I,J) = ZERO
  380    CONTINUE
  400 CONTINUE
C
      IF (OUT) WRITE (IW,9008) SQCDF
      IF (OUT) WRITE (IW,9028) ROTWT
      RETURN
C
 9008 FORMAT(5X,'SQCDF FROM ROTB= ',1P,E15.8)
 9028 FORMAT(5X,'ROTWT USED IN ROTB= ',1P,E15.8)
      END
C*MODULE GVB     *DECK SCHORT
      SUBROUTINE SCHORT(U,NROW,NCOL,NDIM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION U(NDIM,1)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, THRESH=1.0D-10)
C
C     ----- SCHMIDT ORTHONORMALIZE THE EIGENVECTORS IN THE
C           TRANSFORMED REPRESENTATION -----
C
      DO 160 J = 1,NCOL
         SUM = ZERO
         DO 120 K = 1,J
            SUM = DDOT(NROW,U(1,K),1,U(1,J),1)
            IF (K .EQ. J) GO TO 140
            CALL DAXPY(NROW,-SUM,U(1,K),1,U(1,J),1)
  120    CONTINUE
  140    IF (SUM .LT. THRESH) GO TO 180
         SUM = SQRT(SUM)
         CALL DSCAL(NROW,ONE/SUM,U(1,J),1)
  160 CONTINUE
      RETURN
C
  180 IF (MASWRK) THEN
         WRITE (IW,9028) J
         WRITE (IW,9008) SUM
      END IF
      CALL ABRT
      STOP
C
 9008 FORMAT(9H NORM IS ,E13.6)
 9028 FORMAT(1H1/5X,12HEIGENVECTOR ,I4,5X,14HHAS ZERO NORM  /)
      END
C*MODULE GVB     *DECK SELGV1
      SUBROUTINE SELGV1(VECTOR,EIG,L1,IHAM,RSTRCT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL RSTRCT
C
      DIMENSION VECTOR(L1,*),EIG(L1)
      DIMENSION ICH(MXAO)
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
C     ----- SORT THE NEW ORBITALS BASED ON OVERLAP WITH OLD,
C           OR BASED ON LOWEST EIGENVECTORS -----
C
      DO 100 K = 1,L1
         ICH(K) = 0
  100 CONTINUE
C
      IF(.NOT.RSTRCT) THEN
         K = 0
         DO 120 KK = 1,NORB
            K = K+1
            IF (NCONF(KK) .NE. IHAM) GO TO 120
            ICH(K) = 1
  120    CONTINUE
      END IF
C
      CALL SELGV2(VECTOR,EIG,L1,L1,L1,ICH)
      RETURN
      END
C*MODULE GVB     *DECK SELGV2
      SUBROUTINE SELGV2(VECTOR,EIG,NROW,NCOL,NDIM,ICH)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION VECTOR(NDIM,NCOL),EIG(NCOL),ICH(NCOL)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
C     ----- RESTRICT ORBITAL ROTATIONS -----
C     SORT THE ROTATION MATRIX SO THAT THE MAXIMUM COMPONENTS
C     ARE LOCATED ON THE DIAGONAL.  PERMUTE THE EIGENVALUES
C     ACCORDINGLY.  ADJUST THE PHASE OF THE TRANSFORMATION AS WELL.
C
      DO 120 J = 1,NCOL
         KMAX = J
         RMAX = ABS(VECTOR(J,J))
         IF (ICH(J) .NE. 0) RMAX = -EIG(J)
         DO 100 K = J,NCOL
            TEM = ABS(VECTOR(J,K))
            IF (ICH(J) .NE. 0) TEM = -EIG(K)
            IF (TEM .LE. RMAX) GO TO 100
            RMAX = TEM
            KMAX = K
  100    CONTINUE
         IF(J.NE.KMAX) THEN
            TEM       = EIG(J)
            EIG(J)    = EIG(KMAX)
            EIG(KMAX) = TEM
            CALL DSWAP(NROW,VECTOR(1,KMAX),1,VECTOR(1,J),1)
         END IF
         IF(VECTOR(J,J).LT.ZERO) CALL DSCAL(NROW,-ONE,VECTOR(1,J),1)
  120 CONTINUE
      RETURN
      END
C*MODULE GVB     *DECK UPOCBS
      SUBROUTINE UPOCBS(HOLD,TRANS,E,EIG,L1,IHAM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      DIMENSION TRANS(L1,1),HOLD(L1,1),E(1),EIG(1)
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
C     ----- ADDS THE OCBSE CONTRIBUTIONS TO THE NEW COEFFICIENTS -----
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      DO 160 JBF = 1,L1
         FAC = ZERO
         DO 100 K = 1,NORB
            IF (JBF .NE. K) GO TO 100
               FAC = ONE
               IF (NCONF(K)-IHAM) 160,120,160
  100    CONTINUE
  120    CONTINUE
         E(JBF) = EIG(JBF)
         DO 140 I = 1,L1
            HOLD(I,JBF) = HOLD(I,JBF)*FAC +TRANS(I,JBF)
  140    CONTINUE
  160 CONTINUE
      RETURN
      END
C*MODULE GVB     *DECK VIRSET
      SUBROUTINE VIRSET(TRANS,H,H2,D,L1,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      DIMENSION TRANS(L1,*),H(L2),H2(L2),D(L2)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /GVBENG/ EK,VIR,TKIN(25),TPLUSV(25),ASYM
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA NONE/4HNONE/
C
C     ----- ONE ELECTRON ENERGY CONTRIBUTIONS -----
C     -H2- HOLDS H=T+V, -H- HOLDS T INTEGRALS.
C
      CALL DAREAD(IDAF,IODA,H2,L2,11,0)
      IF(MPCTYP.EQ.NONE) CALL DAREAD(IDAF,IODA,H ,L2,13,0)
      DO 100 I = 1,NHAM
         CALL DENSHL(TRANS,D,I,L1,L2)
         IF(MPCTYP.EQ.NONE) TKIN(I) = TRACEP(H,D,L1)
         TPLUSV(I) = TRACEP(H2,D,L1)
  100 CONTINUE
      RETURN
      END
C*MODULE GVB     *DECK ZEROH
      SUBROUTINE ZEROH(H,E,IHAM,IA,L1,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      DIMENSION H(L2),E(L1),IA(L1)
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- ZERO OUT MATRIX ELEMENTS CONNECTING ORBITALS IF THE
C           IHAM'TH FOCK OPERATOR WITH ALL OTHER FOCK OPERATORS -----
C
      DO 200 J = 1,NORB
         IF (NCONF(J) .EQ. IHAM) GO TO 200
         DO 180 K = 1,L1
            IJ = IA(J) + K
            IF(K.GT.J) IJ = IA(K) + J
            H(IJ) = ZERO
  180    CONTINUE
         JJ = IA(J) + J
         H(JJ) = E(J)
  200 CONTINUE
      RETURN
      END
C*MODULE GVB     *DECK BUFDEN
      SUBROUTINE BUFDEN(TRANS,DA,L1,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      COMMON /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
C
      DIMENSION TRANS(L1,1),DA(L2)
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00)
C
C     ----- CALCULATE THE REDUCED DENSITY MATRIX -----
C     ----- FOR THE BUFFER REGION                -----
C
      CALL VCLR(DA,1,L2)
C
      IJ = 0
      DO 280 I = 1,L1
         DO 270 J = 1,I
            IJ = IJ + 1
            SUM = ZERO
            DO 260 K = 1,NBUFMO
               SUM = SUM+F(NCONF(K))*TRANS(I,K)*TRANS(J,K)
  260       CONTINUE
            DA(IJ) = SUM*TWO
  270    CONTINUE
  280 CONTINUE
C
      RETURN
      END
