C 11 APR 08 - MWS - SYNCHRONIZE DFGRID COMMON BLOCK
C  7 DEC 07 - MWS - SCFIN: BONDING ANALYSIS RUN SHOULD FORCE ON VVOS
C 20 AUG 07 - SHY - REDUCE PRINT OUT FOR MD OF QM AND QM/EFP1
C 24 MAR 07 - JAB - SCFIN: USE DIIS IN NON-ABELIAN GROUPS FOR RHF TOO
C 22 DEC 06 - ST,NK,MC,DGF,TN - LC EXCHANGE, FIX FMO/DIIS BUG, FMP/EFP
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 10 JUL 06 - GMC - SCFIN: KILL JOB IF BOTH DIIS AND SOSCF SELECTED
C  8 MAY 06 - MWS - SCFIN: READ VVOS KEYWORD
C 22 FEB 06 - TN  - ENUC: ADD EFPFMO COMMON BLOCK
C 17 JAN 06 - HL  - MODIFY PARALLEL EFP CALCULATIONS
C 14 NOV 05 - DGF - PAD COMMON BLOCKS ENRGYS AND SCFOPT
C 19 SEP 05 - MWS - SCFIN: MP2/CI USE CONV=1D-6 EVEN IF NOT GRADIENT
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 30 APR 05 - DGF - ADD DIRTHR OPTION, FIX GDDI DEADLOCK, SYNCH OPTGRD
C 14 MAR 05 - HL  - ENUC: CALL EXREP INSTEAD OF PAULIR
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 22 FEB 05 - MWS - REMOVE SOME EXTRANEOUS DEBUGGING CODE
C 13 FEB 05 - HMN,JS - ENUC: EFP PARALLELIZATION, EFP MP2 DISPERSION
C  5 FEB 05 - DGF - DMTX,SOGRAD,SOTRAN: RUN SERIAL IF SMALL BASIS SET,
C                   SOGRAD,SOGRRO,SONEWT: TOL DEPENDS ON TARGET CONVHF
C 10 NOV 04 - KI  - DMTX,SOGRAD,SOTRAN: PARALLELIZE SOSCF CONVERGER
C  7 SEP 04 - GDF - DIRFCK: REMOVE POPLE INTEGRAL INDEXING
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C  7 APR 04 - JMR - ENUC: NEEDED GETFM IN EFP2 IF BLOCK
C 16 JAN 04 - MWS - SYNCH COMMON BLOCK DFGRID
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  3 SEP 03 - SPW - MVOS: ADD CIS ARGUMENT TO TWOEI CALL
C  3 JUL 03 - MWS - MVOS: WARN IF ROHF DOESN'T REMOVE ALL SINGLY OCC E-
C 23 MAY 03 - MWS - SCFIN: ABORT GRADIENTS AND MVOQ/ACAVO RUNS
C 12 DEC 02 - MWS - MVOS: SYNCH UP HSTAR ARGUMENTS
C  7 AUG 02 - HT  - USE DYNAMIC MEMORY FOR EFP PAULI REPULSION
C 22 MAY 02 - KRG - DIIS: REDUCE EQUATION SIZE IF SINGULARITY FOUND
C 17 APR 02 - KRG - DIIS: CHANGE PROCEDURE FOR UHF DIIS
C 24 JAN 02 - MWS - SCFIN: READ ACAVO OPTION, CHANGE SOSCF/DIIS LOGICS
C 16 NOV 01 - MWS - MVOS: IMPLEMENT HALF-FILLED VALENCE OPTION
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C 25 JUN 01 - MWS - SCFIN: INPUT FLOATING POINT DENSITY CONVERGENCE
C 13 JUN 01 - MWS - SCFIN: DIIS MADE THE DEFAULT FOR GRID DFT RUNS
C  3 MAY 01 - HL  - ENUC: USE INITIAL GUESS ORBS TO IDENTIFY BUFFER MOS
C 20 FEB 01 - MWS - PAD PAULMO COMMON
C 11 OCT 00 - PB  - CHANGED COMMON BLOCK PCMPAR AND PCMCHG
C 11 JUN 00 - MWS - OVLSEL: FIX FOR SPHERICAL HARMONICS
C 10 APR 00 - MWS - REMOVE STATIC MEMORY FROM COSMO COMMONS
C 25 MAR 00 - KKB/LNB - ENUC: INCLUDE COSMO NUC-SURFACE CHARGE TERM
C 16 FEB 00 - VK  - ENUC: IT IS NECESSARY TO IDENTIFY LMO-S IN QM/MM
C                   OVLSEL: RESTRICT FROZEN MOS IF MFRZ
C  9 APR 99 - MWS - SCFIN: DON'T FORCE DIIS IF USING SOSCF
C 13 MAR 99 - KRG - DIRFCK: ADD CODE TO DO DIRECT DFT
C 13 FEB 99 - KRG - MVOS,MVOORD: CHANGES TO HELP PRESERVE SYMMETRY
C  9 JAN 99 - MWS - DIRFCK: SAVE STATEMENT TO WORK AROUND XLF V5 PROBLEM
C 12 NOV 98 - GDF - SYMH: REMOVE USE OF ISOOUT ROUTINE
C 27 OCT 98 - MAF - LEVELS,OVLSEL,SOGRAD,SOTRAN: ALLOW SPHER.HARM.
C  6 MAY 98 - MWS - DIRFCK: REMOVE DENSITY CRITERION TO CUTOFF TEST
C  6 JAN 98 - MWS - CHANGE CALL TO RACLOS
C 20 NOV 97 - GMC - SOHESS: CAP LARGEST HESSIAN ELEMENT AT HTOL.
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C  2 SEP 97 - PND - ENUC: PRINT 9300 ONLY ON MASTER NODE
C 14 AUG 97 - MWS - ENUC: ADD DECLARATION FOR POLNAM
C 19 JUL 97 - MWS - SONEWT: ABORT CORRECTLY IF RUNNING IN PARALLEL
C 16 JUL 97 - GNM - ENUC: CHANGES FOR FRAGONLY RUNS, PRINT CHANGES
C 18 MAR 97 - MWS - ENUC: ADD PCM CORRECTION TO NUCLEAR REPULSION
C 18 DEC 96 - TLW - DIRFCK: CODE FOR RESPONSE FOCK-LIKE MATRICES
C  9 DEC 96 - JHJ - ENUC: CALL TO PAULIR
C  3 DEC 96 - MWS - MVOS: ALLOW FOR DIRECT AND PARALLEL COMPUTATION
C 17 OCT 96 - GMC - SONEWT: SCALE BIG STEPS, ABORT HUGE ONES
C 18 SEP 96 - MWS - ENUC: PASS DYNAMIC MEMORY TO EREPUL
C 13 JUN 96 - VAG - DELETE THE OLD 'SCF' ROUTINE
C 24 MAY 96 - MWS - INVOKE ENUCQ AS A SUBROUTINE CALL
C 12 MAY 96 - MWS - SONEWT: FIX BROADCAST MISTAKE
C 18 APR 96 - GMC - SONEWT: SKIP IF GRADIENT IS VERY SMALL
C  5 MAR 96 - MWS - CHANGE SHELL SYMMETRY PACKING COMMONS
C 14 SEP 95 - WC  - ENUC: ADDITIONAL FRAGMENT OUTPUT ADDED
C 18 AUG 95 - MWS - SCFIN: TURN SOSCF OFF FOR VERY SMALL BASES
C 11 AUG 95 - XL  - ENUC: ADD SCREENING TERM
C 26 JUL 95 - GMC - IMPLEMENT SOSCF FOR UHF,ROHF,GVB,
C                   NEW CODE SOGRRO,GVBGRD,GVBHES,GVBTRN ADDED
C 27 MAR 95 - MWS - SCFIN: MAKE SO GRADIENT TOLERANCE AN INPUT VALUE,
C                   CORRECT F,ALPHA,BETA FOR ROHF WITH NO OPEN SHELLS
C  8 MAR 95 - GMC - SOGRAD,SOHESS,SONEWT,SOTRAN ADDED FOR SOSCF
C  5 MAR 95 - MWS - DIRFCK: CHANGE SCREENING OF INDIVIDUAL INTS
C  1 FEB 95 - MWS - ENUC,EXTRAP: CHANGES FOR MOROKUMA ANALYSIS
C 15 DEC 94 - TLW - ENUC: ELECTRIC FIELD CONTRIBUTION ADDED
C 12 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 30 JUN 94 - JHJ - SCF: CALL LMOX FOR LMO-MP2
C 27 MAY 94 - SPW - ENUC: FRAGMENT CORRECTIONS NOW DONE IN EFPCMP
C 27 MAY 94 - PND - ENUC: FRAGMENT ENERGY CORRECTIONS ADDED
C 18 APR 94 - MWS - SCFIN: DIIS DEFAULT OFF FOR SMALL BASES
C 12 APR 94 - MWS - MVOS: USE JACOBI DIAGONALIZATION
C  5 APR 94 - SPW - SCF: REMOVE INDUCTION OF SCRF
C 25 MAR 94 - MWS - SCF: RESTORE ROHF CANONICAL MO-S AFTER A ROMP2
C 10 MAR 94 - FRJ - SCFIN: ALLOW READ OF INITIAL EXPONENT SCALE FACTOR
C 10 AUG 93 - MWS - SCFIN: TEST THE NO ARRAY FOR SENSIBILITY
C 16 JUL 93 - MWS - INCREASE MAXIMUM CI ROOTS TO 100
C  1 APR 93 - PND - ENUC: PRINT NUCLEAR CHARGE TERM
C 11 MAR 93 - FJ  - INCLUSION OF VIRIAL ANALYSIS
C  8 MAR 93 - JHJ - ENUC: REMOVE CONVERSION FACTOR
C 28 JAN 93 - MR  - DIISER: THRESHOLD TESTS IN ERROR COMPUTATION
C 13 JUN 92 - MWS - DIIS INTERPOLATION FORCED TO RUN SEQUENTIALLY,
C                   NEW ROUTINE DIISER FOR ERROR MATRIX CALCULATION
C  2 APR 92 - MWS,TLW - COMMON ENRGYS MADE PURE FLOATING POINT
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C  2 MAR 92 - MWS - DELETE LVLSHF, ADD LEVELS AND OVLSEL
C 19 FEB 92 - JHJ - SCFIN: DIIS NOT DEFAULT, AND NO DEM FOR MPC WFN.
C 11 FEB 92 - JHJ - ENUC: ADD MOPAC NRE VIA COMMON ENUCLR.
C  7 JAN 92 - TLW - MAKE WRITES PARALLEL; ADD COMMON PAR
C 17 DEC 91 - MWS - DIIS: REMOVE CONVERGENCE TEST
C 13 NOV 91 - JHJ - SCF:ENERGS(1) IS FILLED BEFORE RHFCL CALL.
C 31 OCT 91 - MWS - SHLDEN: FIX FOR DECSTN'S F77 3.0 COMPILER
C 31 OCT 91 - JHJ - SCF: WIDENED LOOP SO THAT COMMON ENRGYS IS FILLED.
C 23 OCT 91 - JHJ - INTRODUCED DIELEC AND ICALC TO COMMON /ZRFPAR/.
C 16 OCT 91 - JHJ - ENUC:DEFINED GRAD WHICH IS USED INSTEAD OF OPTIMZ.
C 11 OCT 91 - MWS - DIIS: ERROR MATRIX PASSED IN, FOCK NOT PUT ON DAF,
C                   DIRFCK,SHLDEN: NEW ROUTINES ADDED FOR DIRECT SCF,
C                   SCFIN: READS DIRECT SCF KEYWORDS.
C 18 AUG 91 - MWS - DIIS: GVB CASE READS ERROR IN AO BASIS FROM DAF,
C                   SCFIN: READ DEM INPUT, ELIMINATE BTHRSH INPUT,
C                   MAKFV,MAKEIJ: CLEANUP
C 15 AUG 91 - TLW - TRMAT: EXPAND FOR F AND G FUNCTIONS
C 14 AUG 91 - TLW - SYMH: AND RHR: EXPANDED FOR F AND G FUNCTIONS
C 23 JUL 91 - JHJ - ENUC:ADDED CALL TO DNUCPO.
C 17 JUL 91 - JHJ - ENUC:ADDED DRG'S ZERNER FIELD CHANGES.
C  2 JUL 91 - JHJ - ENUC:ADDED DRG'S MULT.POLE CHANGES TO GMS ENUC.
C                   MOVED ENUCD,ENUCQ,ENUCO TO EFINT.
C 22 MAY 91 - MWS - SCFIN: FIX COUNTING OF OPEN SHELL ELECTRONS
C  5 FEB 91 - MWS - SCFIN: FIX GVB COUPLINGS FOR NO CORE ORBITAL CASE
C 12 SEP 90 - MWS - INTRODUCE MXATM
C 15 AUG 90 - DRG - MOVE ENUC CALL TO SCF FROM THE SEPARATE ROUTINES.
C 14 AUG 90 - DRG - EFP REACTION FIELD LOOP ADDED TO SCF.
C  7 AUG 90 - TLW - ADD CF AND CG TO COMMON NSHEL
C 19 JUN 90 - MWS - SCFIN: TEST $SCF ELECTRON COUNT WITH INPUT CHARGE
C  9 JUN 90 - STE - EXTRAP,LVLSHF: MOVE SAVE STATEMENTS
C 16 DEC 89 - MWS - SCFIN: SET KCORB CORRECTLY
C 23 OCT 89 - MWS - DELETE /FUNCT/ SAVE
C 10 JUL 89 - MWS - SCFIN: FILL F,ALPHA,BETA,ETC FOR RHF CASE.
C  9 MAY 89 - MWS - SCFIN: FILL F,ALPHA,BETA FOR HIGH SPIN ROHF,
C                   CORRECT F,ALPHA,BETA FOR TRIPLET GVB CASE (THE OLD
C                   PARAMETERS WERE WRONG IF DEGENERATE OPEN SHELLS),
C                   CHANGE DEFAULT CICOEF FROM 0.9,0.2 TO 0.9,-0.2!
C 12 FEB 89 - MWS - DON'T CALL WFNMP2 IF SCF DOESN'T CONVERGE.
C 27 JAN 89 - MWS - CALL WFNMP2 FROM ROUTINE SCF
C 15 DEC 88 - MWS - RENAME /SPD/ TO /SYMSPD/, DELETE FROM SCF AND SYMH
C 14 NOV 88 - MWS - ADD ROUTINE FOR MVOS
C 10 AUG 88 - MWS - MXSH,MXGSH,MXGTOT FROM 120,10,440 TO 1000,30,5000
C 19 JUL 88 - MWS - DEFAULT BTHRSH FROM 3 TO 10, BMAX AT LEAST 0.25
C 29 JUN 88 - JAM - FIX DIIS FOR GVB-PP(1) CASE
C  8 MAY 88 - MWS - UP AO-S FROM 255 TO 2047
C 22 APR 88 - MWS - USE GVBWFN AND SCFWFN COMMONS IN SCFIN.
C 17 APR 88 - MWS - CLEAN UP DIIS, CHANGES IN SCF AND SCFIN RELATED
C                   TO DOING ROHF IN THE UHF CODE.
C  1 APR 88 - MWS - CALL NEW ROUTINE SLVLEQ TO SOLVE DIIS EQNS.
C  4 FEB 88 - MWS - ADD MORE COUPLING CASES TO ROUTINE SCFIN
C 14 NOV 87 - STE - MOVE SCFIN FROM INPUT; SAVE /FUNCT/
C 27 APR 87 - STE - DMTX: REMOVE IA ARRAY
C 24 NOV 86 - MWS - CHANGE RSHIFT TO RRSHFT
C 19 NOV 86 - STE - DIIS: CHANGE IER TO IERR; PARAMETERS, NO NAV USE
C 12 AUG 86 - MWS - BRING DIIS ROUTINES FROM MODULE DIIGAM,
C                   USE LINPACK INSTEAD OF IMSL,  REMOVE BRENDA
C                   LAM'S DIIGAM MODULE OF APRIL 1986.
C  7 JUL 86 - JAB - RENAME ROUTINE SHIFTR TO LVLSHF
C  7 JUL 86 - MWS - SANITIZE FLOATING POINT CONSTANTS
C 23 SEP 85 - STE - EXTRAP,SHIFTR: SAVE VARIABLES. GENERIC FUNCT.S
C                   DAMPD: INSERT GO TO 360 BEFORE 340 CONTINUE
C 17 FEB 85 - MWS - ADD THE ROUTINE CLENMO
C 28 JUN 84 - STE - EXPAND RECORD 2 (ENERGY VALUES) OF DAF 10 IN SCF
C 18 JAN 84 - STE - CHANGE LENGTH OF /FUNCT/ IN SCF
C 27 DEC 83 - STE - REPLACE SCPRD WITH DDOT IN EXTRAP; WFN TO GAMESS
C  5 NOV 83 - STE - DELETE CALL TO BUFINT IN SCF
C 19 APR 83 - MWS - ICALP,ICBET IN EXTRAP
C 22 MAR 83 - MWS - SPLIT SCFA INTO SCFLIB AND RHFUHF
C                   THESE ARE THE UTILITY ROUTNIES USED FOR HF METHODS
C  9 MAR 83 - MWS - CONVERGE IF DELTA-E IS TEENY FOR RHF,UHF
C  4 MAR 83 - MWS - PUNCH ORBITALS UNLESS NPUNCH IS 0
C  4 OCT 82 - MWS - CONVERT TO THE IBM
C 20 JUN 82 - MWS - NPRINT=8 TO GET FULL PRINTOUT IN RHFCL,UHFOP,
C                   GVBOUT,RHFOP - EIGENVALUES, DENSITY SQUELCHED
C
C*MODULE SCFLIB  *DECK ACAVO
      SUBROUTINE ACAVO(PACAVO,VEC,EIG,DCORE,DVAL,
     *                 JCORE,JVAL,KCORE,KVAL,
     *                 Q,SCR,IWRK,BUFFX,BUFFI,
     *                 NOPK,NINTMX,L0,L1,L2,L3,NFROZ)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DOUBLE PRECISION JCORE(L2),JVAL(L2),KCORE(L2),KVAL(L2)
      DIMENSION PACAVO(6)
      DIMENSION DCORE(L2),DVAL(L2),VEC(L1,L1),EIG(L1),
     *          Q(L1,L1),SCR(L1,8),IWRK(L1),
     *          BUFFX(NINTMX),BUFFI(NINTMX)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00, SMALL=1.0D-10)
C
C     ----- CONSTRUCT GENERAL CORRELATION ADAPTED VIRTUAL ORBITALS -----
C     PROGRAMMED BY MIKE SCHMIDT AT IOWA STATE IN DECEMBER 2001.
C     THE CONCEPT IS FROM LAIMUTIS BYTAUTAS AND KLAUS RUEDENBERG.
C
C     THE ACAVO'S DIAGONALIZE THE VIRTUAL SPACE OF THE OPERATOR
C       F(CORR) =  A (KINETIC ENERGY)       + B (NUCLEAR ATTRACTION)
C                + C (COUL.REP.BY CORE ELS) + D (COUL.REP.BY VAL.ELS)
C                + E (EXCH.REP.BY CORE ELS) + F (EXCH.REP.BY VAL.ELS)
C     THE CASE A=B=C=D=E=0, AND F=-1 CORRESPONDS TO WHITTEN ORBITALS,
C         J.L.WHITTEN  J.CHEM.PHYS.  56, 5458-5466(1972)
C
C     THE FIRST -NFROZ- ORBITALS IN VEC AND EIG WILL BE UNALTERED.
C
      NCORE = NUMCOR()
      NVAL  = NFROZ - NCORE
      NVIRT = L0 - NFROZ
      LVIRT = NFROZ+1
      IF (MASWRK) WRITE(IW,9000) NVIRT,NCORE,NVAL,PACAVO
C
C     ----- REORDER VIRTUALS BY SYMMETRY TYPE -----
C     THIS LETS THE JACOBI DIAGONALIZATION MORE EASILY PRESERVE SYMMETRY
C
      CALL DAREAD(IDAF,IODA,VEC,L3,15,0)
      CALL DAREAD(IDAF,IODA,EIG,L1,17,0)
      CALL MVOORD(VEC,EIG,SCR(1,1),SCR(1,2),SCR(1,3),SCR(1,4),
     *            Q,DCORE,L0,L1,L2,L3,NFROZ)
C
C         GET CORE AND VALENCE DENSITY MATRICES, BY SUBTRACTING
C         THE CORE DENSITY FROM THE TOTAL DENSITY OF ANY TYPE SCF.
C
      CALL DENDD1(DVAL,DCORE,L2)
C
      DO I=1,NCORE
         SCR(I,1) = TWO
      ENDDO
      CALL DMTX(DCORE,VEC,SCR,NCORE,L1,L1)
C
      CALL VSUB(DCORE,1,DVAL,1,DVAL,1,L2)
C
C         MAKE CORE AND VALENCE COULOMB AND EXCHANGE OPERATORS.
C         THIS IS CODED FOR INTEGRALS ON DISK ONLY, AT PRESENT.
C
      IF(NOPK.NE.1) THEN
         IF(MASWRK) WRITE(IW,9010)
         CALL ABRT
      END IF
      CALL HSTARJK(DCORE,DVAL,JCORE,JVAL,KCORE,KVAL,
     *             BUFFX,BUFFI,NINTMX,IA,NOPK,L1,L2)
C
C         CONSTRUCT THE OPERATOR IN THE AO BASIS IN -DCORE-,
C         USING MATRIX -DVAL- AS SCRATCH STORAGE.
C
      FT  = PACAVO(1)
      FV  = PACAVO(2)
      FJC = PACAVO(3)
      FJV = PACAVO(4)
      FKC = PACAVO(5)
      FKV = PACAVO(6)
C
C         READ -H- AND -T- TO OBTAIN -V- BY SUBTRACTION.
C
      CALL DAREAD(IDAF,IODA,DVAL ,L2,11,0)
      CALL DAREAD(IDAF,IODA,DCORE,L2,13,0)
      CALL VSUB(DCORE,1,DVAL,1,DVAL,1,L2)
C
      CALL DSCAL(L2,FT ,        DCORE,1)
      CALL DAXPY(L2,FV ,DVAL ,1,DCORE,1)
      CALL DAXPY(L2,FJC,JCORE,1,DCORE,1)
      CALL DAXPY(L2,FJV,JVAL ,1,DCORE,1)
      CALL DAXPY(L2,FKC,KCORE,1,DCORE,1)
      CALL DAXPY(L2,FKV,KVAL ,1,DCORE,1)
C
C     ----- TRANSFORM TO VIRTUAL MO BASIS, AND MAKE ACAVO'S -----
C     USE JACOBI DIAGONALIZATION TO TRY TO AVOID MIXING DEGENERATE MO-S
C
      CALL TFTRI(DVAL,DCORE,VEC(1,LVIRT),SCR,NVIRT,L1,L1)
      DO 300 I=1,(NVIRT*NVIRT+NVIRT)/2
        IF (ABS(DVAL(I)) .LT. SMALL) DVAL(I)=ZERO
  300 CONTINUE
      IGERR=0
      KDSAVE = KDIAG
      KDIAG = 3
      CALL GLDIAG(L1,NVIRT,NVIRT,DVAL,SCR,EIG(LVIRT),
     *            VEC(1,LVIRT),IGERR,IWRK)
      KDIAG = KDSAVE
      IF(IGERR.GT.0) CALL ABRT
      CALL DAREAD(IDAF,IODA,Q,L3,15,0)
      CALL TFSQB(VEC(1,LVIRT),Q(1,LVIRT),Q,NVIRT,L1,L1)
C
      CALL DAWRIT(IDAF,IODA,VEC,L3,15,0)
      CALL DAWRIT(IDAF,IODA,EIG,L1,17,0)
      RETURN
C
 9000 FORMAT(/1X,'GENERATING',I5,' GENERAL CORRELATION ADAPTED',
     *           ' VIRTUAL ORBITALS...'/
     *        1X,'THERE ARE',I5,' FROZEN CORE ORBITALS, AND',I5,
     *           ' FROZEN VALENCE ORBITALS'/
     *        1X,'ACAVO PARAMETERS=',6F10.5)
 9010 FORMAT(1X,'CANNOT FORM J,K MATRICES USING SUPERMATRIX INTEGRALS'/
     *       1X,'PLEASE RESUBMIT WITH NOPK=1')
      END
C
C*MODULE MP2     *DECK CLENMO
C     ***********************************************
      SUBROUTINE CLENMO(VEC,NAO,NMO,TOLZ,TOLE,IW,OUT)
C     ***********************************************
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL OUT
C
      DIMENSION VEC(NAO,NMO)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
C      GENERATE SYMMETRY ADAPTED MO-S BY SETTING MO COEFFICIENTS
C      SMALLER THAN TOLZ TO ZERO, AND SETTING MO COEFFICIENTS
C      WHOSE MAGNITUDE AGREES TO WITHIN TOLE EXACTLY EQUAL.
C      AUTHOR - MIKE SCHMIDT, FEBRUARY 6, 1985
C
      NZER=0
      NEQU=0
      DO 400 J=1,NMO
         DO 300 I=1,NAO
            VAL1=ABS(VEC(I,J))
            IF(VAL1.EQ.ZERO) GO TO 300
C                              ZERO THE TEENY ONES
            IF(VAL1.GT.TOLZ) GO TO 100
               VEC(I,J)=ZERO
               NZER=NZER+1
               GO TO 300
  100       CONTINUE
C                              EQUAL MAGNITUDES FOR THE NONZERO
            IF(I.EQ.NAO) GO TO 300
            FAC=ONE/VAL1
            KMIN=I+1
            DO 200 K=KMIN,NAO
               VAL2=ABS(VEC(K,J))
               TEST=ABS(FAC*(VAL1-VAL2))
               IF(TEST.GT.TOLE) GO TO 200
               IF(TEST.EQ.ZERO) GO TO 200
                  VEC(K,J) = SIGN(VAL1,VEC(K,J))
                  NEQU=NEQU+1
  200       CONTINUE
  300    CONTINUE
  400 CONTINUE
      IF(OUT) WRITE(IW,900) NZER,NEQU
      RETURN
C
  900 FORMAT(1X,I5,' MO COEFS WERE ZEROED,',I5,' WERE EQUATED')
      END
C*MODULE SCFLIB  *DECK DAMPD
      SUBROUTINE DAMPD(DE,DEP,DEAVG,DAMP,ACURCY,DIFF,DIFFP,DMPTLC)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00, PT5=0.5D+00, TWO=2.0D+00, FOUR=4.0D+00)
      PARAMETER (PT2=0.2D+00, PT25=0.25D+00, FAC=16.0D+00)
C
      DAMPO = DAMP
      ETEST = ACURCY*ACURCY
      IF (ABS(DE) .LT. ETEST .AND. ABS(DEP) .LT. ETEST) GO TO 300
      IF (ABS(DE) .LT. ETEST) GO TO 320
      IF (ABS(DEP) .LT. ETEST) GO TO 340
      IF ((DIFFP-DIFF) .LT. ZERO) GO TO 100
      IF ( ABS(DE) .GE. ACURCY .OR. DE .GT. ZERO) GO TO 100
C
C     ----- CONVERGED -----
C
      DAMP = DAMP/FAC
      GO TO 280
C
  100 CONTINUE
      IF ( DE .GT. ZERO) GO TO 200
      IF (DEP .GT. ZERO) GO TO 180
      IF ( DE .GT. DEP) GO TO 140
C
C     ----- DE < 0. , DEP < 0. , DE < DEP -----
C
      IF ( ABS(DE) .LT. TWO*DEAVG) GO TO 120
      DAMP = FAC*MAX(DAMP,DEAVG)
      GO TO 280
C
  120 IF ( ABS(DE) .GT. PT5*DEAVG) GO TO 280
      DAMP = DAMP/FAC
      GO TO 280
C
  140 CONTINUE
C
C     ----- DE < 0. , DEP < 0. , DE > DEP -----
C
      IF (DE .GT. PT25*DEP) GO TO 160
      DAMP = (DE/DEP)**2*MAX(DAMP,DEAVG)
      GO TO 280
C
  160 DAMP = DAMP/FAC
      GO TO 280
C
  180 CONTINUE
C
C     ----- DE < 0. , DEP > 0. -----
C
      DAMP = FOUR*MAX(DAMP,DEAVG)
      IF (-DE .GT. DEAVG) DAMP = DAMP*FAC
      IF (-DE+DEP .GE. DEAVG) GO TO 280
      DAMP = DAMP/FAC
      GO TO 280
C
  200 CONTINUE
      IF (DEP .GT. ZERO) GO TO 220
C
C     ----- DE > 0. , DEP < 0. -----
C
      DAMP = FOUR*MAX(DAMP,DEAVG)
      IF (DE .GT. PT5*DEAVG) DAMP = DAMP*FAC
      IF (DE-DEP .GE. PT2*DEAVG) GO TO 280
      DAMP = DAMP/FAC
      GO TO 280
C
  220 CONTINUE
C
C     ----- DE > 0. , DEP > 0. -----
C
      DAMP = FOUR*MAX(DAMP,DEAVG)
      IF (DE .LT. FOUR*DEP) GO TO 240
      DAMP = FAC*MAX(DAMP,DEAVG)
      GO TO 280
C
  240 IF (DE .GT. PT25*DEP) GO TO 260
      DAMP = DAMP/FAC
      GO TO 280
C
  260 DAMP = (DE/DEP)**2*MAX(DAMP,DEAVG)
  280 CONTINUE
C
C     ----- IF THE DENSITY CONVERGENCE WORSENED - MAKE SURE
C           THAT THE DAMPING CAN'T DECREASE -----
C
      IF ((DIFFP-DIFF) .LT. ZERO) DAMP = MAX(DAMP,DAMPO)
      GO TO 360
C
  300 CONTINUE
C
C        DE < ETEST AND DEP < ETEST
      DAMP = DAMP/FAC
      GO TO 360
C
  320 CONTINUE
C        DE < ETEST  DEP > ETEST
      DAMP = DAMP/FAC
      GO TO 360
C
  340 CONTINUE
C        DEP < ETEST  DE > ETEST
      DAMP = DAMPO
      IF (DE .GT. ZERO) DAMP = MAX(TWO*DAMP,DMPTLC)
C
  360 CONTINUE
      RETURN
      END
C*MODULE SCFLIB  *DECK DDIFF
      SUBROUTINE DDIFF(D0,D1,L2,DIFF)
C
C     ----- CALCULATE THE LARGEST ABSOLUTE CHANGE IN THE
C           DENSITY MATRIX -----
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION D0(2),D1(2)
C
      DIFF = 0.0D+00
      DO 100 I = 1,L2
      DIFIJ = ABS(D0(I)-D1(I))
      IF (DIFIJ .GT. DIFF) DIFF = DIFIJ
  100 CONTINUE
      RETURN
      END
C*MODULE SCFLIB  *DECK DIIS
      SUBROUTINE DIIS(SCFTYP,IW,ITDIIS,Q,FCKA,FCKB,ERR,ERRB,WRK,
     *                A,X,IPVT,B,IODIIS,SCR,L1,L2,L3,MAXIT,MAXIT2,
     *                MAXIO,ERDIIS,NOTOPN,MASOUT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL NOTOPN,DBG,GOPARR,DSKWRK,MASWRK,GPSAVE,MASOUT
C
      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
      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
C
      DIMENSION Q(L3),FCKA(L2),FCKB(L2),ERR(L3),WRK(L3),IODIIS(MAXIO),
     *          A(MAXDII,MAXDII),X(MAXIT),IPVT(MAXIT),B(MAXIT2),SCR(L1),
     *          ERRB(L3)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      DATA DBUG,DBUGME/8HDEBUG   ,8HDIIS    /
      DATA UHF/8HUHF     /
C
C                 THIS ROUTINE PERFORMS THE
C        *** DIRECT INVERSION IN THE ITERATIVE SUBSPACE ***
C        INTERPOLATION USING THE CURRENT AND PREVIOUS FOCK
C        MATRICES TO SPEED UP THE CONVERGENCE OF THE SCF.
C
C        THE METHOD IS P.PULAY, J.COMPUT.CHEM., 3, 556-560 (1982).
C        BUT SEE ALSO P.PULAY, CHEM.PHYS.LETT., 73, 393-398 (1980),
C        T.P.HAMILTON, P.PULAY, J.CHEM.PHYS., 84, 5728-5734 (1986).
C
C        WRITTEN BY BRENDA LAM IN APRIL 1986 AT U. HOUSTON
C        SIMPLIFIED BY MWS IN SEPTEMBER 1991 AT NDSU.
C
C     ON ENTRY, -FCKA- AND -FCKB- ARE THE CURRENT FOCK MATRICES.
C               -ERR- AND -ERRB- ARE THE CURRENT ERROR MATRICES,
C               EXPRESSED IN THE AO BASIS.
C     ON EXIT, -ERR- AND -ERRB- ARE DESTROYED, WHILE -FCKA- AND -FCKB-
C               WILL CONTAIN INTERPOLATED FOCK MATRICES (IF ITDIIS.GT.2)
C
C     ITDIIS.EQ.1 MEANS DIIS IS NOT YET TURNED ON
C     ITDIIS.EQ.2 MEANS DIIS HAS BEEN INITIATED
C     ITDIIS.GT.2 MEANS A DIIS INTERPOLATION OF ITDIIS-1
C                 FOCK MATRICES HAS BEEN PERFORMED.
C
      DBG = (NPRINT.EQ.5  .OR.  EXETYP.EQ.DBUG  .OR.  EXETYP.EQ.DBUGME)
     *      .AND. MASWRK
C
C     DIIS DATA IS SAVED ON DIRECT FILE -IRAF- IN CYCLIC ORDER,
C                 1,2,3,...,MAXDII = LAST MAXDII ERROR MATRICES
C        MAXDII + 1,2,3,...,MAXDII = LAST MAXDII ERROR MATRICES (BETA)
C      2*MAXDII + 1,2,3,...,MAXDII = LAST MAXDII ALPHA FOCK MATRICES
C      3*MAXDII + 1,2,3,...,MAXDII = LAST MAXDII BETA  FOCK MATRICES
C
      LRAFEA= 0
      LRAFEB= MAXDII
      LRAFFA= MAXDII*2
      LRAFFB= MAXDII*3
C
C     ----- IF WE ARE ON THE LAST SCF CYCLE, SHUT DIIS DOWN -----
C
      IF(ITDIIS.GE.MAXIT) THEN
         IF(.NOT.NOTOPN) CALL RACLOS(IRAF,'DELETE')
         NOTOPN = .TRUE.
         RETURN
      END IF
C
      IF(NOTOPN) THEN
         NOTOPN = .FALSE.
         IRAF   = 20
         NDAF20 = MAXIO
         LEN20  = L3
         CALL RAOPEN(IRAF,IODIIS,0,NDAF20,LEN20,-5)
      END IF
C
C     ----- PUT ERROR MATRIX INTO CONSISTENT O.N.B. -----
C     PULAY USES S**-1/2, BUT HERE WE USE Q, Q OBEYS Q-DAGGER*S*Q=I
C     E-ORTH = Q-DAGGER * E * Q, FCKA IS USED AS A SCRATCH -L1- VECTOR
C
      CALL DCOPY(L3,ERR,1,WRK,1)
      CALL TFSQU(ERR,WRK,Q,SCR,L1,L1)
      IF(SCFTYP.EQ.UHF)THEN
        CALL DCOPY(L3,ERRB,1,WRK,1)
        CALL TFSQU(ERRB,WRK,Q,SCR,L1,L1)
      END IF
      IF(DBG) THEN
         WRITE(IW,*) 'ERROR MATRIX IN O.N.B.'
         CALL PRSQ(ERR,L1,L1,L1)
      END IF
C
C     ----- START DIIS PROCEDURE IF ERDIIS < ETHRSH -----
C
      IMAX   = IDAMAX(L3,ERR,1)
      ERDIIS = ABS(ERR(IMAX))
      IF(SCFTYP .EQ. UHF) THEN
        JMAX   = IDAMAX(L3,ERRB,1)
        ERDIIS = ERDIIS + ABS(ERRB(JMAX))
      END IF
      IF(ERDIIS.GT.ETHRSH  .AND.  ITDIIS.EQ.1) RETURN
C
C     ----- SAVE THE CURRENT FOCK MATRIX -----
C
      LFCKA = LRAFFA + MOD(ITDIIS-1,MAXDII) + 1
      LFCKB = LRAFFB + MOD(ITDIIS-1,MAXDII) + 1
      CALL RAWRIT(IRAF,IODIIS,FCKA,L2,LFCKA,0)
      IF(SCFTYP.EQ.UHF) CALL RAWRIT(IRAF,IODIIS,FCKB,L2,LFCKB,0)
C
C     ----- SAVE THE CURRENT ERROR MATRIX -----
C
      LERR = LRAFEA + MOD(ITDIIS-1,MAXDII) + 1
      CALL RAWRIT(IRAF,IODIIS,ERR,L3,LERR,0)
      IF(SCFTYP.EQ.UHF)THEN
        LERR = LRAFEB + MOD(ITDIIS-1,MAXDII) + 1
        CALL RAWRIT(IRAF,IODIIS,ERRB,L3,LERR,0)
      END IF
C
      ITDIIS = ITDIIS+1
C
C        WHEN RUNNING IN PARALLEL, THE O(N**2) COMMUNICATIONS
C        COSTS IN -RAREAD- BELOW ARE THE SAME MAGNITUDE AS
C        THE COMPUTATION IN -DDOT- AND -DAXPY-.
C        IT SEEMS WISEST JUST TO LET THE MASTER DO THIS I/O,
C        THEN BROADCAST THE INTERPOLATION RESULTS AT THE END.
C
      GPSAVE = GOPARR
      GOPARR = .FALSE.
C
C        FORM LAST ROW (ROW ITDIIS) OF TRIANGULAR MATRIX B,
C        FIRST DOING THE DIAGONAL AND THE FIRST ELEMENT OF THE ROW.
C        THE TRANSPOSE IN PULAY'S STEP (2) IS JUST A SIGN CHANGE
C        BECAUSE ERROR MATRICES ARE ANTISYMMETRIC BY CONSTRUCTION.
C        THE SIGN IS IGNORED BECAUSE IT DOESN'T AFFECT THE LINEAR
C        EQUATION'S SOLUTION.
C
      BJJ=DDOT(L3,ERR,1,ERR,1)
      IF(SCFTYP .EQ. UHF) THEN
        BJJ = (BJJ + DDOT(L3,ERRB,1,ERRB,1))/2.0D+00
      END IF
C
      IF(ITDIIS.EQ.2) THEN
         IF (MASOUT) WRITE(IW,900)
         CALL VCLR(B,1,MAXIT2)
         B(1)=ZERO
         B(2)=-ONE
         B(3)= BJJ
         GOPARR = GPSAVE
         RETURN
      ELSE
         J1   = (ITDIIS*ITDIIS-ITDIIS)/2 + 1
         JJ   = (ITDIIS*ITDIIS+ITDIIS)/2
         B(J1)= -ONE
         B(JJ)= BJJ
      END IF
C
C        THE REST OF THE BIJ'S.
C        NOTE THAT WE ONLY COMPUTE THE NA-2 VALUES NEXT TO THE
C        DIAGONAL, SO THAT THE OTHER ELEMENTS (EXCEPT THE FIRST
C        COLUMN) ARE RANDOM VALUES.
C
      NA=MIN(ITDIIS,MAXDII)
      NAM1=NA-1
      IJ=JJ-NA+1
      DO 120 IX=2,NAM1
         LERR = LRAFEA + MOD(ITDIIS-NA+IX-2,MAXDII) + 1
         CALL RAREAD(IRAF,IODIIS,WRK,L3,LERR,0)
         BIJ=DDOT(L3,ERR,1,WRK,1)
         IF(SCFTYP.EQ.UHF)THEN
           LERR = LRAFEB + MOD(ITDIIS-NA+IX-2,MAXDII) + 1
           CALL RAREAD(IRAF,IODIIS,WRK,L3,LERR,0)
           BIJ=(BIJ+DDOT(L3,ERRB,1,WRK,1))/2.0D+00
         END IF
         IJ=IJ+1
         B(IJ)=BIJ
  120 CONTINUE
C
      IF(DBG) THEN
         WRITE(IW,*) 'DIIS EQUATIONS'
         CALL PRTRI(B,ITDIIS)
      END IF
C
C     ----- SET UP AND SOLVE THE SET OF DIIS LINEAR EQUATIONS -----
C
  180 CONTINUE
      CALL DIISEQ(A,NA,B,JJ,ITDIIS)
      CALL VCLR(X,1,NA)
      X(1)=-ONE
      IERR=0
      CALL SLVLEQ(A,X,IPVT,NA,NA,0,IERR)
      IF(IERR.NE.0) THEN
         NA = NA - 1
         IF(MASWRK) WRITE(IW,*)
     *      'REDUCING DIIS EQUATION SIZE BY 1 FOR NUMERICAL STABILITY'
         IF (NA .NE. 1) GO TO 180
         IF(MASWRK) WRITE(IW,*) 'DIIS EQUATIONS ARE SINGULAR, BOMBING'
         CALL ABRT
      END IF
C
      IF(DBG) THEN
         WRITE(IW,*) 'DIIS INTERPOLATION COEFFICIENTS'
         WRITE(IW,910) (X(I),I=1,NA)
      END IF
C
C     ----- FORM THE INTERPOLATED ALPHA FOCK MATRIX -----
C     THE INTERPOLATION IS SUMMATION C(I)*F(I)
C
      CALL VCLR(FCKA,1,L2)
      DO 220 IX=2,NA
         LFCKA = LRAFFA + MOD(ITDIIS-NA+IX-2,MAXDII) + 1
         CALL RAREAD(IRAF,IODIIS,WRK,L2,LFCKA,0)
         CI = X(IX)
         CALL DAXPY(L2,CI,WRK,1,FCKA,1)
  220 CONTINUE
      IF(DBG) THEN
         WRITE(IW,*) 'DIIS INTERPOLATED ALPHA FOCK MATRIX IS'
         CALL PRTRIL(FCKA,L1)
      END IF
C
C     ----- FORM THE INTERPOLATED -UHF- BETA FOCK MATRIX -----
C
      IF(SCFTYP.EQ.UHF) THEN
         CALL VCLR(FCKB,1,L2)
         DO 240 IX=2,NA
            LFCKB = LRAFFB + MOD(ITDIIS-NA+IX-2,MAXDII) + 1
            CALL RAREAD(IRAF,IODIIS,WRK,L2,LFCKB,0)
            CI = X(IX)
            CALL DAXPY(L2,CI,WRK,1,FCKB,1)
  240    CONTINUE
         IF(DBG) THEN
            WRITE(IW,*) 'DIIS INTERPOLATED BETA FOCK MATRIX IS'
            CALL PRTRIL(FCKB,L1)
         END IF
      END IF
C
      GOPARR=GPSAVE
      IF(GOPARR) THEN
         CALL DDI_BCAST(1005,'F',FCKA,L2,MASTER)
         IF(SCFTYP.EQ.UHF)
     *   CALL DDI_BCAST(1006,'F',FCKB,L2,MASTER)
      END IF
      RETURN
C
  900 FORMAT(10X,'* * *   INITIATING DIIS PROCEDURE   * * *')
  910 FORMAT(1X,5F13.8)
      END
C*MODULE SCFLIB  *DECK DIISEQ
      SUBROUTINE DIISEQ(A,NA,B,NTT,NB)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(NA,NA), B(NTT)
C
C        FILL SQUARE MATRIX A WITH SYMMETRIC MATRIX B SKIPPING
C        ROWS (AND COLUMNS) BEGINNING FROM SECOND ROW (COLUMN).
C        WRITTEN BY BRENDA LAM - APRIL 1986
C
      NSKP1=NB-NA+1
C
      IB=0
      JA=0
      DO 120 J=1,NB
         IA=0
         DO 110 I=1,J
            IB=IB+1
            IF(I.EQ.1 .OR. I.GT.NSKP1) THEN
               IF(J.EQ.1 .OR. J.GT.NSKP1) THEN
                  IA=IA+1
                  IF(IA.EQ.1) JA=JA+1
                  BIB=B(IB)
                  A(IA,JA)=BIB
                  A(JA,IA)=BIB
               END IF
            END IF
  110    CONTINUE
  120 CONTINUE
      RETURN
      END
C*MODULE SCFLIB  *DECK DIISER
      SUBROUTINE DIISER(F,D,S,ERR,WRK,L1,L2,IPASS,NPASS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION F(L1,L1),D(L2),S(L1,L1),ERR(L1,L1),WRK(L1)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT,PARR
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXSEQ=150, MXROWS=5)
      PARAMETER (ZERO=0.0D+00)
C
C     ----- CALCULATE THE DIIS ERROR MATRIX, ERR=FDS-SDF -----
C     THIS ROUTINE IS GENERAL FOR RHF, ROHF, AND UHF.
C
C     SINCE EACH OF F,D,S ARE SYMMETRIC, THIS MEANS CALCULATE FDS,
C     AND THEN SUBTRACT THE TRANSPOSE FROM THAT RESULT.
C     BEFORE ENTRY, -F- AND -S- MUST BE EXPANDED TO SQUARE STORAGE.
C     -IPASS- EQUAL 1 OR 2 MEANS ALPHA VERSUS BETA ERROR TERMS.
C
C     ----- INITIALIZATION FOR PARALLEL WORK -----
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT  = -1
      L2CNT = -1
      PARR = GOPARR  .AND.  L1.GT.MXSEQ
C
      DO 310 I = 1,L1,MXROWS
         IIMAX = MIN(L1,I+MXROWS-1)
C
C     ----- GO PARALLEL! -----
C
         IF(PARR) THEN
            IF (NXT) THEN
               L2CNT = L2CNT + 1
               IF (L2CNT.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
               IF (NEXT.NE.L2CNT) THEN
                  IF(IPASS.EQ.1) THEN
                     DO 010 II=I,IIMAX
                        CALL VCLR(ERR(II,1),L1,L1)
  010                CONTINUE
                  END IF
                  GO TO 310
               END IF
            ELSE
               IPCOUNT = IPCOUNT + 1
               IF (MOD(IPCOUNT,NPROC).NE.0) THEN
                  IF(IPASS.EQ.1) THEN
                     DO 020 II=I,IIMAX
                        CALL VCLR(ERR(II,1),L1,L1)
  020                CONTINUE
                  END IF
                  GO TO 310
               END IF
            END IF
         END IF
C
C             FIRST CALCULATE -F- TIMES -D- ONE ROW AT A TIME.
C             SINCE -F- IS SYMMETRIC, WE ACTUALLY USE F-DAGGER.
C
         DO 300 II=I,IIMAX
            KL = 0
            DO 140 L = 1,L1
               LM1 = L-1
               DUM = ZERO
               FDUM = F(L,II)
               IF (LM1.GT.0) THEN
                  DO 100 K = 1,LM1
                     KL = KL+1
                     WRK(K) = WRK(K)+D(KL)*FDUM
                     DUM = DUM+D(KL)*F(K,II)
  100             CONTINUE
               END IF
               KL = KL+1
               WRK(L) = DUM+D(KL)*FDUM
  140       CONTINUE
C
C             THEN TAKE THAT ROW OF -FD- TIMES EVERY COLUMN IN -S-.
C             THE UHF BETA ERROR IS ACCUMULATED TO THE ALPHA ERROR.
C
            IF(IPASS.EQ.1) THEN
               DO 200 J = 1,L1
                  DUM = DDOT(L1,WRK,1,S(1,J),1)
                  IF(ABS(DUM).LT.1.0D-15) DUM=ZERO
                  ERR(II,J) = DUM
  200          CONTINUE
            ELSE
               DO 210 J = 1,L1
                  DUM = DDOT(L1,WRK,1,S(1,J),1)
                  IF(ABS(DUM).LT.1.0D-15) DUM=ZERO
                  ERR(II,J) = ERR(II,J) + DUM
  210          CONTINUE
            END IF
  300    CONTINUE
  310 CONTINUE
C
      IF(PARR) THEN
         IF(IPASS.EQ.NPASS) CALL DDI_GSUMF(1008,ERR,L1*L1)
         IF(NXT) CALL DDI_DLBRESET
      END IF
C
C        FINISH BY SUBTRACTING THE TRANSPOSE OF F*D*S FROM ITSELF
C
      IF(IPASS.EQ.NPASS) CALL SUBMT(ERR,L1)
C
      RETURN
      END
C*MODULE SCFLIB  *DECK DIRFCK
      SUBROUTINE DIRFCK(SCFTYP,IA,DA,FA,DB,FB,GHONDO,L2,NINT,NXYZ)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL IANDJ,KANDL,SAME,OUT,UROHF,LCFLAG,LRINT
C
      DIMENSION IA(*),DA(*),FA(*),DB(*),FB(*),GHONDO(*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXAO=8192)
C
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /ERIOUT/ ISH,JSH,KSH,LSH,LSTRI,LSTRJ,LSTRK,LSTRL
      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,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MISC  / IANDJ,KANDL,SAME
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      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 /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJX,IJ,KL,IJKL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
C
      DATA HALF /0.5D+00/
      DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF     /
      DATA GVB/8HGVB     /
C
C-NEXT STATEMENT IS FOR VARIOUS IBM XLF 3.X AND 5.X COMPILERS-
C
      SAVE IJN,KLN
C
C     ----- FORM FOCK OPERATOR DIRECTLY FROM INTEGRALS -----
C     THIS ROUTINE WAS PIECED TOGETHER FROM QOUT AND HSTAR
C     BY FRANK JENSEN AT ODENSE UNIVERSITY IN MARCH 1990.
C     SCF FUNCTIONS BESIDES RHF ADDED BY MWS IN AUGUST 1991.
C
C     NOTE THAT OFF-DIAGONAL ELEMENTS WILL NEED TO BE HALVED LATER.
C
      HFSCAL=DFTTYP(3)
      CSCALT=1.0D+00
      IF(LCFLAG) THEN
       IF(LRINT) THEN
       HFSCAL=1.0D+00
       CSCALT=0.0D+00
       ELSE
       HFSCAL=0.0D+00
       CSCALT=1.0D+00
       ENDIF
      ENDIF
C
      UROHF = SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
      CUTINT = CUTOFF
C
      IF(SCFTYP.EQ.GVB) THEN
         NSHL = NHAM
         IF(NCO.GT.0) NSHL=NSHL-1
      END IF
C
      SAME  = ISH .EQ. KSH .AND. JSH .EQ. LSH
      IANDJ = ISH .EQ. JSH
      KANDL = KSH .EQ. LSH
C
      MINI = KMIN(ISH)
      MINJ = KMIN(JSH)
      MINK = KMIN(KSH)
      MINL = KMIN(LSH)
      MAXI = KMAX(ISH)
      MAXJ = KMAX(JSH)
      MAXK = KMAX(KSH)
      MAXL = KMAX(LSH)
      LOCI = KLOC(ISH)-MINI
      LOCJ = KLOC(JSH)-MINJ
      LOCK = KLOC(KSH)-MINK
      LOCL = KLOC(LSH)-MINL
C
      IJN = 0
      JMAX = MAXJ
      DO 360 I = MINI,MAXI
         I_INDEX = (I-MINI)*LSTRI + 1
         IF (IANDJ) JMAX = I
         DO 340 J = MINJ,JMAX
            IJ_INDEX = (J-MINJ)*LSTRJ + I_INDEX
            IJN = IJN+1
            LMAX = MAXL
            KLN = 0
            DO 320 K =  MINK,MAXK
               IJK_INDEX = (K-MINK)*LSTRK + IJ_INDEX
               IF (KANDL) LMAX = K
               DO 300 L = MINL,LMAX
                  KLN = KLN+1
                  IF(SAME .AND. KLN.GT.IJN) GO TO 340
                  IJKL_INDEX = (L-MINL)*LSTRL + IJK_INDEX
C
                  VAL = GHONDO( IJKL_INDEX )
                  IF(ABS(VAL).LT.CUTINT) GO TO 300
                  NINT = NINT + 1
C
                  II = LOCI+I
                  JJ = LOCJ+J
                  KK = LOCK+K
                  LL = LOCL+L
C
                  IF (II .GE. JJ) GO TO 100
                  N = II
                  II = JJ
                  JJ = N
  100             IF (KK .GE. LL) GO TO 120
                  N = KK
                  KK = LL
                  LL = N
  120             IF (II-KK) 140,160,180
  140             N = II
                  II = KK
                  KK = N
                  N = JJ
                  JJ = LL
                  LL = N
                  GO TO 180
  160             IF (JJ .LT. LL) GO TO 140
  180             CONTINUE
C
                  IF(II.EQ.JJ) VAL = VAL*HALF
                  IF(KK.EQ.LL) VAL = VAL*HALF
                  IF(II.EQ.KK  .AND.  JJ.EQ.LL) VAL = VAL*HALF
C
C      WE NOW HAVE EVERYTHING READY, PUT INTEGRALS INTO FOCK MATRIX
C
                  NIJ = IA(II)+JJ
                  NKL = IA(KK)+LL
                  NIK = IA(II)+KK
                  NIL = IA(II)+LL
                  NJK = IA(JJ)+KK
                  NJL = IA(JJ)+LL
                  IF(JJ.LT.KK) NJK = IA(KK)+JJ
                  IF(JJ.LT.LL) NJL = IA(LL)+JJ
C
                  VAL2 = VAL+VAL
                  VAL4 = VAL2+VAL2
C
C       NXYZ DISTINGUISHES CODE FOR RHF OR RHF RESPONSE EQUATIONS
C
C       HFSCAL DISTINGUISHES CODE FOR HF OR DFT
C
                  IF(HFSCAL .EQ. 1.0D+00.AND..NOT.LCFLAG) THEN
                  IF(SCFTYP.EQ.RHF) THEN
                     IF(NXYZ.EQ.1) THEN
                        FA(NIJ) = FA(NIJ)+VAL4*DA(NKL)
                        FA(NKL) = FA(NKL)+VAL4*DA(NIJ)
                        FA(NIK) = FA(NIK)- VAL*DA(NJL)
                        FA(NIL) = FA(NIL)- VAL*DA(NJK)
                        FA(NJK) = FA(NJK)- VAL*DA(NIL)
                        FA(NJL) = FA(NJL)- VAL*DA(NIK)
                     ELSE
                        NIJ = (NIJ-1)*NXYZ+1
                        NKL = (NKL-1)*NXYZ+1
                        NIK = (NIK-1)*NXYZ+1
                        NIL = (NIL-1)*NXYZ+1
                        NJK = (NJK-1)*NXYZ+1
                        NJL = (NJL-1)*NXYZ+1
                        DO 210 IXYZ=0,NXYZ-1
                          FA(NIJ+IXYZ) = FA(NIJ+IXYZ)+VAL4*DA(NKL+IXYZ)
                          FA(NKL+IXYZ) = FA(NKL+IXYZ)+VAL4*DA(NIJ+IXYZ)
                          FA(NIK+IXYZ) = FA(NIK+IXYZ)- VAL*DA(NJL+IXYZ)
                          FA(NIL+IXYZ) = FA(NIL+IXYZ)- VAL*DA(NJK+IXYZ)
                          FA(NJK+IXYZ) = FA(NJK+IXYZ)- VAL*DA(NIL+IXYZ)
                          FA(NJL+IXYZ) = FA(NJL+IXYZ)- VAL*DA(NIK+IXYZ)
  210                   CONTINUE
                     END IF
                  END IF
C
                  IF(UROHF) THEN
                     DUM = VAL4*(DA(NKL)+DB(NKL))
                     FA(NIJ) = FA(NIJ)+DUM
                     FB(NIJ) = FB(NIJ)+DUM
                     DUM = VAL4*(DA(NIJ)+DB(NIJ))
                     FA(NKL) = FA(NKL)+DUM
                     FB(NKL) = FB(NKL)+DUM
                     FA(NIK) = FA(NIK)-VAL2*DA(NJL)
                     FB(NIK) = FB(NIK)-VAL2*DB(NJL)
                     FA(NIL) = FA(NIL)-VAL2*DA(NJK)
                     FB(NIL) = FB(NIL)-VAL2*DB(NJK)
                     FA(NJK) = FA(NJK)-VAL2*DA(NIL)
                     FB(NJK) = FB(NJK)-VAL2*DB(NIL)
                     FA(NJL) = FA(NJL)-VAL2*DA(NIK)
                     FB(NJL) = FB(NJL)-VAL2*DB(NIK)
                  END IF
C
                  IF(SCFTYP.EQ.GVB) THEN
                     IF(NCO.GT.0) THEN
                        FA(NIJ) = FA(NIJ)+VAL4*DA(NKL)
                        FA(NKL) = FA(NKL)+VAL4*DA(NIJ)
                        FA(NIK) = FA(NIK)-VAL *DA(NJL)
                        FA(NIL) = FA(NIL)-VAL *DA(NJK)
                        FA(NJK) = FA(NJK)-VAL *DA(NIL)
                        FA(NJL) = FA(NJL)-VAL *DA(NIK)
                        IOFF1 = L2
                        IOFF2 = L2
                     ELSE
                        IOFF1 = 0
                        IOFF2 = 0
                     END IF
                     DO 220 IFO = 1,NSHL
                        FA(NIJ+IOFF1)    = FA(NIJ+IOFF1)
     *                              + VAL4*DA(NKL+IOFF2)
                        FA(NKL+IOFF1)    = FA(NKL+IOFF1)
     *                              + VAL4*DA(NIJ+IOFF2)
                        FA(NIK+IOFF1+L2) = FA(NIK+IOFF1+L2)
     *                              + VAL *DA(NJL+IOFF2)
                        FA(NIL+IOFF1+L2) = FA(NIL+IOFF1+L2)
     *                              + VAL *DA(NJK+IOFF2)
                        FA(NJK+IOFF1+L2) = FA(NJK+IOFF1+L2)
     *                              + VAL *DA(NIL+IOFF2)
                        FA(NJL+IOFF1+L2) = FA(NJL+IOFF1+L2)
     *                              + VAL *DA(NIK+IOFF2)
                        IOFF1 = IOFF1+L2+L2
                        IOFF2 = IOFF2+L2
  220                CONTINUE
                  END IF
                  ELSE
C
C                 DFT/LC-DFT CASE
C
                  IF(SCFTYP.EQ.RHF) THEN
                     IF(NXYZ.EQ.1) THEN
                        FA(NIJ) = FA(NIJ)+VAL4*DA(NKL)*CSCALT
                        FA(NKL) = FA(NKL)+VAL4*DA(NIJ)*CSCALT
                        FA(NIK) = FA(NIK)- VAL*DA(NJL)*HFSCAL
                        FA(NIL) = FA(NIL)- VAL*DA(NJK)*HFSCAL
                        FA(NJK) = FA(NJK)- VAL*DA(NIL)*HFSCAL
                        FA(NJL) = FA(NJL)- VAL*DA(NIK)*HFSCAL
                     ELSE
                        NIJ = (NIJ-1)*NXYZ+1
                        NKL = (NKL-1)*NXYZ+1
                        NIK = (NIK-1)*NXYZ+1
                        NIL = (NIL-1)*NXYZ+1
                        NJK = (NJK-1)*NXYZ+1
                        NJL = (NJL-1)*NXYZ+1
                        DO 310 IXYZ=0,NXYZ-1
                          FA(NIJ+IXYZ) = FA(NIJ+IXYZ)+VAL4*DA(NKL+IXYZ)
     *                                                          *CSCALT
                          FA(NKL+IXYZ) = FA(NKL+IXYZ)+VAL4*DA(NIJ+IXYZ)
     *                                                          *CSCALT
                          FA(NIK+IXYZ) = FA(NIK+IXYZ)- VAL*DA(NJL+IXYZ)
     *                                                          *HFSCAL
                          FA(NIL+IXYZ) = FA(NIL+IXYZ)- VAL*DA(NJK+IXYZ)
     *                                                          *HFSCAL
                          FA(NJK+IXYZ) = FA(NJK+IXYZ)- VAL*DA(NIL+IXYZ)
     *                                                          *HFSCAL
                          FA(NJL+IXYZ) = FA(NJL+IXYZ)- VAL*DA(NIK+IXYZ)
     *                                                          *HFSCAL
  310                   CONTINUE
                     END IF
                  END IF
C
                  IF(UROHF) THEN
                     DUM = VAL4*(DA(NKL)+DB(NKL))*CSCALT
                     FA(NIJ) = FA(NIJ)+DUM
                     FB(NIJ) = FB(NIJ)+DUM
                     DUM = VAL4*(DA(NIJ)+DB(NIJ))*CSCALT
                     FA(NKL) = FA(NKL)+DUM
                     FB(NKL) = FB(NKL)+DUM
                     FA(NIK) = FA(NIK)-VAL2*DA(NJL)*HFSCAL
                     FB(NIK) = FB(NIK)-VAL2*DB(NJL)*HFSCAL
                     FA(NIL) = FA(NIL)-VAL2*DA(NJK)*HFSCAL
                     FB(NIL) = FB(NIL)-VAL2*DB(NJK)*HFSCAL
                     FA(NJK) = FA(NJK)-VAL2*DA(NIL)*HFSCAL
                     FB(NJK) = FB(NJK)-VAL2*DB(NIL)*HFSCAL
                     FA(NJL) = FA(NJL)-VAL2*DA(NIK)*HFSCAL
                     FB(NJL) = FB(NJL)-VAL2*DB(NIK)*HFSCAL
                  END IF
C
                  IF(SCFTYP.EQ.GVB) THEN
                    WRITE(IW,*)'DIRECT GVB DFT NOT IMPLEMENTED'
                    CALL ABRT()
C                   I COULDN'T FIGURE OUT THE DO 220 LOOP  :-)
                  END IF
                  END IF
C
  300          CONTINUE
  320       CONTINUE
  340    CONTINUE
  360 CONTINUE
      RETURN
      END
C*MODULE SCFLIB  *DECK DMTX
      SUBROUTINE DMTX(D,V,X,M,N,NDIM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION D(*),V(NDIM,M),X(M)
C
      PARAMETER (ZERO=0.0D+00, MXSEQ=150)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,PARR
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARR = GOPARR  .AND.  N.GT.MXSEQ
      DO 120 I = 1,N
C
C     --- GO PARALLEL! ---
C
         IF(PARR) THEN
            IF(MOD(I,NPROC).NE.ME) GO TO 120
         END IF
C
         IJ = I*(I-1)/2
         DO 110 J = 1,I
            IJ = IJ + 1
            DUM = ZERO
            DO 100 K = 1,M
               DUM = DUM+X(K)*V(I,K)*V(J,K)
  100       CONTINUE
            D(IJ) = DUM
  110    CONTINUE
  120 CONTINUE
C
C     --- GO PARALLEL! ---
C     COLLECT ROWS OF -D- FROM EACH PROCESS, THEN BROADCAST TO ALL
C
      IF(PARR) THEN
         DO 130 I=1,N
            IFROM = MOD(I,NPROC)
            IF(IFROM.EQ.MASTER) GO TO 130
            IJ=(I*I-I)/2+1
            IF(MASWRK) THEN
               CALL DDI_RECV(D(IJ),8*I,IFROM)
            ELSE
               IF(IFROM.EQ.ME) CALL DDI_SEND(D(IJ),8*I,MASTER)
            END IF
  130    CONTINUE
         CALL DDI_BCAST(1011,'F',D,(N*N+N)/2,MASTER)
      END IF
C
      RETURN
      END
C*MODULE SCFLIB  *DECK ENUC
      DOUBLE PRECISION FUNCTION ENUC(N,Z,C)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,EFLDL,SCREEN,LEFNUC,IFGRAD
      LOGICAL LPRINT
C
      DOUBLE PRECISION MOROKM, MD
C
      PARAMETER (MXFRG=1050, MXDFG=5, MXPT=2000,
     *           MXFGPT=12000, MXATM=2000, MXAO=8192,
     *           MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION Z(N),C(3,N)
      DIMENSION ELTOT(10)
C
      COMMON /EFLDC / EVEC(3),EFLDL
      COMMON /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      COMMON /EFPFMO/ IEFPFMO,ISWENUC,REPNUCEFP,EFPEFP,RUNEFP
      COMMON /ENUCLR/ ENUCLR
      COMMON /FMCOM / XX(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 /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OPTGRD/ Q(3*MXATM),ES,QQ(20),
     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGH
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB(MXFRG),
     *                NPBF(MXFRG),NTMO
      COMMON /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FE,UNZOR,PORT
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IRPPCM,IEF,IP_F,NFMOPCM
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCINP / VLAMB,SCREEN
      COMMON /SIMDAT/ NACC,NREJ,IGOMIN,NRPA,IBWM,NACCT,NREJT,NRPAT,
     *                NPRTGO,IDPUNC,IGOFLG
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /ZMTALT/ NZMAT2,NZVAR2,NVAR2,NZMTRD,ICOORD
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
C
C
C     COSMO INFORMATION
C
      PARAMETER (NPPA=1082)
      LOGICAL ISEPS,USEPS
      COMMON /ISEPS / ISEPS,USEPS
      COMMON /COSMO1/ SE2,SECORR,ETOTS,CDUM,QVCOSMO,
     *                CSPOT(NPPA),ICORR,ITRIPO,ITRIP2,ITRIP3,ITRIP4,
     *                NATCOS,NPS,ITERC
      COMMON /COSMO3/ COSZAN(NPPA),CORZAN(3,NPPA)
C
      PARAMETER (ZERO=0.0D+00)
C
      DATA MOROKM /8HMOROKUMA/, MD /8HMD      /
      DATA NONE/4HNONE/
C
C     --- MOPAC NUCLEAR REPULSION ENERGY IS CALCULATED IN MHCORE  ---
C
      IF(MPCTYP.NE.NONE) THEN
         REPNUC=ENUCLR
         GO TO 800
      END IF
C
C     --- MOROKUMA ANALYSIS ---
C
      IF (RUNTYP.EQ.MOROKM) THEN
         NN=N
         CALL STENUC(NN,NST,NED)
      ELSE
         NST = 2
         NED = N
      END IF
C
      REPNUC = ZERO
      DUM2 = ZERO
      DUM3 = ZERO
      DUM4 = ZERO
      DUM5 = ZERO
      ENEL=ZERO
      IF (N.EQ.1) GO TO 310
C
C     CHECKING IF IT IS AN SCREEN CALCULATION. IF SO REPULSION TERM
C     IS ZA*ZB*EXP(-LAMBDA*RAB)/RAB
C
      IF (SCREEN) THEN
         DO 301 I = NST,NED
            NI = I-1
            DO 201 J = NST-1,NI
               RR = ZERO
               DO 101 K = 1,3
                   RR = RR+(C(K,I)-C(K,J))**2
  101          CONTINUE
               REPNUC = REPNUC+Z(I)*Z(J)*EXP(-VLAMB*SQRT(RR))/SQRT(RR)
  201       CONTINUE
  301    CONTINUE
      ELSE
C
         DO 300 I = NST,NED
            NI = I-1
            DO 200 J = NST-1,NI
               RR = ZERO
               DO 100 K = 1,3
                  RR = RR+(C(K,I)-C(K,J))**2
  100          CONTINUE
               IF(RR.NE.ZERO) REPNUC = REPNUC+Z(I)*Z(J)/SQRT(RR)
  200       CONTINUE
  300    CONTINUE
      END IF
C
C     NUCLEAR - SAS CHARGE INTERACTION FROM COSMO SOLVENT MODEL
C
      CDUM = ZERO
      IF(ISEPS  .AND.  MPCTYP.EQ.NONE  .AND.  ICORR.EQ.1) THEN
         WRITE(IW,*)'INSIDE SCFLIB, EXTRA NUCLEAR LOOP, N,NPS=',N,NPS
         DO 122 I=1,NPS
            DO 124 J=1,N
               RR=ZERO
               DO 126 K=1,3
                  RR=RR+(CORZAN(K,I)-C(K,J))**2
  126          CONTINUE
               CDUM = CDUM + COSZAN(I)*Z(J)/SQRT(RR)
  124       CONTINUE
  122    CONTINUE
         WRITE(IW,*)'NUCLEAR CHARGE (REPNUC), NO NUCLEAR-SAS:',REPNUC
         WRITE(IW,*)'NUCLEAR-SAS CONTRIBUTION (CDUM):',CDUM
         REPNUC=REPNUC+CDUM
         WRITE(IW,*)'TOTAL VALUE OF NUCLEAR-CHARGE COMPONENT:',REPNUC
      END IF
C
C           ADD NUCLEAR CONTRIBUTION FROM ELECTRIC FIELD
C
  310 CONTINUE
      IF (EFLDL) THEN
         ANUCF = ZERO
         DO 330 J = 1,N
            DO 320 I = 1,3
               ANUCF = ANUCF - EVEC(I)*C(I,J)*Z(J)
  320       CONTINUE
  330    CONTINUE
         REPNUC = REPNUC + ANUCF
      END IF
C
C           ADD NUCLEAR CONTRIBUTION FROM PCM MODEL
C
      IF(IPCM.EQ.1) THEN
         IF(MASWRK) WRITE(IW,9200) UNZ
         REPNUC = REPNUC+UNZ
      END IF
C
C     ---- VARIOUS EFP ADDITIVE CONSTANTS INCLUDED IN THE V-NN ---
C
      ISUM = IEFC+IEFD+IEFQ+IEFO+IEFP+NFRG
      IF(ISUM.EQ.0) GO TO 800
C
      LEFNUC = NFRG.GT.1.AND.ISWENUC.EQ.1
C
      IF (ICOORD .EQ. 4) THEN
         IDUM = NSERCH
         NSERCH = 0
         CALL EFSP
         NSERCH = IDUM
         GO TO 800
      END IF
C
      DO I=1,10
        ELTOT(I) = ZERO
      ENDDO
C
C   CALL TO MOIDBM IDENTIFIES ATOMS WHICH ARE IN THE BUFFER.
C   NOTE THAT THIS USES THE INITIAL GUESS ORBITALS.
C
      IF(NBUFMO.GT.0) CALL MOIDBM(265)
C
C    ATTENTION! WE CALL HERE FFELEC ROUTINE WHICH DOES ENERGY
C               (AND GRADIENTS) FOR ALL EFP-EFP ELECTROSTATIC TERMS
      IF(LEFNUC) THEN
         MAXDER = 0
         IFGRAD = .FALSE.
         CALL DERCHK(MAXDER)
         IF (MAXDER.GT.0) IFGRAD = .TRUE.
         CALL FFELEC(ESNRG,ELTOT(1),ELTOT(2),ELTOT(4),ELTOT(7),
     *               ELTOT(3),ELTOT(5),ELTOT(6),IFGRAD,.FALSE.)
      END IF
C
      IF(IEFC.EQ.1) DUM2 = ENUCC(N,Z,C)
      IF(IEFD.EQ.1) DUM3 = ENUCD(N,Z,C)
      IF(IEFQ.EQ.1) CALL ENUCQ(DUM4,N,Z,C)
      IF(IEFO.EQ.1) DUM5 = ENUCO(N,Z,C)
C
C        NOTE THAT THERE ARE TWO DIFFERENT REPULSION FORMULA FOR (8)
C
      IF(IREP.EQ.1.AND.LEFNUC) THEN
         MXRPTS=0
         DO 400 IFRG=1,NFRG
            MXRPTS = MAX(NRPTS(IFRG),MXRPTS)
  400    CONTINUE
         MXRPTS=MXRPTS+1
         NDFRG2 = (NDFRG*NDFRG+NDFRG)/2
         LEN = MXRPTS*MXRPTS*NDFRG2
         CALL VALFM(LOADFM)
         LAREL = LOADFM + 1
         LCREL = LAREL  + LEN
         LAST  = LCREL  + LEN
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
         CALL DAREAD(IDAF,IODA,XX(LAREL),2*LEN,75,0)
         CALL EREPUL(ELTOT(8),XX(LAREL),XX(LCREL),MXRPTS,NDFRG2)
         CALL RETFM(NEED)
      END IF
C
      IF (NTMO.GT.0 .AND. LEFNUC) THEN
         MXMO=0
         MXBF=0
         DO I = 1, NFRG
            MXMO=MAX(MXMO,NORB(I))
            MXBF=MAX(MXBF,NPBF(I))
         END DO
         MXMO2=(MXMO*MXMO+MXMO)/2
C
         CALL VALFM(LOADFM)
         LPROVEC = LOADFM  + 1
         LFOCKMA = LPROVEC + MXBF*NTMO
         LSMAT   = LFOCKMA + MXMO2*NFRG
         LTMAT   = LSMAT   + MXBF*MXBF
         LWRK    = LTMAT   + MXBF*MXBF
         LSIJ    = LWRK    + MXBF
         LTIJ    = LSIJ    + MXMO*MXMO
         LFASQ   = LTIJ    + MXMO*MXMO
         LFBSQ   = LFASQ   + MXMO*MXMO
         LAST    = LFBSQ   + MXMO*MXMO
         NEED=LAST-LOADFM-1
         CALL GETFM(NEED)
C
         LENPV=MXBF*NTMO
         LENFM=MXMO2*NFRG
         CALL DAREAD(IDAF,IODA,XX(LPROVEC),LENPV,267,0)
         CALL DAREAD(IDAF,IODA,XX(LFOCKMA),LENFM,268,0)
         CALL EXREP(ELTOT(8),XX(LPROVEC),XX(LFOCKMA),XX(LSMAT),
     *              XX(LTMAT),XX(LWRK),XX(LSIJ),XX(LTIJ),
     *              XX(LFASQ),XX(LFBSQ),MXBF,MXMO)
         CALL RETFM(NEED)
      END IF
C
      IF(IEFP2.EQ.1.AND.LEFNUC) THEN
         MXDPTS = 0
         DO 450 IFRG = 1,NFRG
            MXDPTS = MAX(NDPTS(IFRG),MXDPTS)
  450    CONTINUE
         MXDPTS = MXDPTS+1
         NDFRG2 = (NDFRG*NDFRG+NDFRG)/2
         LEN = MXDPTS*MXDPTS*NDFRG2
C
         CALL VALFM(LOADFM)
         LAREL = LOADFM + 1
         LDREL = LAREL + LEN
         LCREL = LDREL + LEN
         LAST  = LCREL + LEN
         LBREL = LAST  + LEN
         LEREL = LBREL + LEN
         LAST1 = LEREL + LEN
         NEED  = LAST1 - LOADFM - 1
         CALL GETFM(NEED)
         CALL DAREAD(IDAF,IODA,XX(LAREL),6*LEN,81,0)
C
         CN=6.0D+00
         CALL EDISPR(ELTOT(9),XX(LAREL),XX(LDREL),XX(LCREL),
     *               MXDPTS,NDFRG2,CN)
         CN=8.0D+00
         CALL EDISPR(ELTOT(10),XX(LAST),XX(LBREL),XX(LEREL),
     *               MXDPTS,NDFRG2,CN)
         CALL RETFM(NEED)
      END IF
C
C     - ONLY THE ELECTROSTATIC TERMS NEED TO SUM UP
C
      DO I=1,10
         ENEL = ENEL + ELTOT(I)
      ENDDO
      REPNUC0 = REPNUC
C
      IF (ISWENUC.EQ.1.AND.IEFPFMO.EQ.1) THEN
        REPNUC = REPNUC+DUM2+DUM3+DUM4+DUM5
        REPNUCEFP = ENEL
      ELSE
        REPNUC = REPNUC+DUM2+DUM3+DUM4+DUM5+ENEL
      END IF
C
      LPRINT=MASWRK  .AND.  NPRTGO.NE.2 .AND. ISWENUC.EQ.1
      IF(RUNTYP.EQ.MD) LPRINT=.FALSE.
      IF(LPRINT) THEN
         WRITE(IW,9008) REPNUC0
         WRITE(IW,9018) DUM2
         WRITE(IW,9028) DUM3
         WRITE(IW,9038) DUM4
         WRITE(IW,9048) DUM5
         IF(LEFNUC) THEN
            WRITE(IW,9118) ELTOT(1)
            WRITE(IW,9119) ELTOT(2)
            WRITE(IW,9121) ELTOT(4)
            WRITE(IW,9124) ELTOT(7)
            WRITE(IW,9120) ELTOT(3)
            WRITE(IW,9122) ELTOT(5)
            WRITE(IW,9123) ELTOT(6)
            WRITE(IW,9125) ELTOT(8)
            WRITE(IW,9127) ELTOT(9)
            WRITE(IW,9128) ELTOT(10)
            WRITE(IW,9095) ENEL
         END IF
         WRITE(IW,9199) REPNUC
      END IF
C
  800 CONTINUE
      ENUC = REPNUC
      RETURN
C
 9008 FORMAT(/1X,'EFP ADDITIVE CONTRIBUTIONS TO THE "NUCLEAR ENERGY":'/
     *       1X,'  -- NUCLEAR CHARGE - NUCLEAR CHARGE   =',F20.10)
 9018 FORMAT(1X,'  -- NUCLEAR CHARGE - EFP CHARGE       =',F20.10)
 9028 FORMAT(1X,'  -- NUCLEAR CHARGE - EFP DIPOLE       =',F20.10)
 9038 FORMAT(1X,'  -- NUCLEAR CHARGE - EFP QUADRUPOLE   =',F20.10)
 9048 FORMAT(1X,'  -- NUCLEAR CHARGE - EFP OCTUPOLE     =',F20.10)
 9095 FORMAT(1X,'  -- TOTAL INTERFRAGMENT INTERACTION   =',F20.10)
 9118 FORMAT(1X,'  -- EFP CHARGE     - EFP CHARGE       =',F20.10)
 9119 FORMAT(1X,'  -- EFP CHARGE     - EFP DIPOLE       =',F20.10)
 9121 FORMAT(1X,'  -- EFP CHARGE     - EFP QUADRUPOLE   =',F20.10)
 9124 FORMAT(1X,'  -- EFP CHARGE     - EFP OCTUPOLE     =',F20.10)
 9120 FORMAT(1X,'  -- EFP DIPOLE     - EFP DIPOLE       =',F20.10)
 9122 FORMAT(1X,'  -- EFP DIPOLE     - EFP QUADRUPOLE   =',F20.10)
 9123 FORMAT(1X,'  -- EFP QUADRUPOLE - EFP QUADRUPOLE   =',F20.10)
 9125 FORMAT(1X,'  -- EFP-EFP EXCHANGE REPULSION ENERGY =',F20.10)
 9127 FORMAT(1X,'  -- EFP-EFP DISPERSION ENERGY (C6)    =',F20.10)
 9128 FORMAT(1X,'  -- EFP-EFP DISPERSION ENERGY (C8)    =',F20.10)
 9199 FORMAT(1X,'THE TOTAL NUCLEAR CHARGE/EFP ENERGY    =',F20.10)
 9200 FORMAT(1X,'NUCLEI-NUCLEAR CHARGE POLARIZATION ADJUSTMENT TO VNN=',
     *        F20.10)
      END
C*MODULE SCFLIB  *DECK EXTRAP
      SUBROUTINE EXTRAP(DE,DAMP,DAMP0,H0,H1,H2,H3,L1,L2,LL2,NDAF,ITERV,
     *                  NCALL,ITYP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION H0(*),H1(*),H2(*),H3(*)
C
      LOGICAL EXTRAH,DAMPH,VSHIFT,EXTPRE,DAMPRE
C
      SAVE DAMPRE,EXTPRE
C
      COMMON /ACONV / RRSHFT,EXTTOL,DMPTOL,VSHTOL,IEXTIN
      COMMON /CONV  / ACURCY,EN,ETOT,EHF,EHF0,DIFF,ITER,ICALP,ICBET
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00)
      PARAMETER (PT5=0.5D+00, TM17=1.0D-17)
      PARAMETER (TOL=1.0D-07, TOL1=1.9D+00, TOL2=0.99D+00)
      PARAMETER (SHRNKF=100.0D+00, DMPMIN=0.01D+00, RSMIN=0.8D+00)
C
      DATA DAMPRE,EXTPRE/.FALSE.,.FALSE./
C
C     ----- DECODE CONVERGENCE METHOD -----
C     USE EITHER DAVIDSON'S DAMPING AND EXTRAPOLATION
C            OR POPLE'S EXTRAPOLATION
C
      EXTRAH = MOD(MCONV,2) .EQ. 1
      DAMPH  = MOD(MCONV,4) .GE. 2
      VSHIFT = MOD(MCONV,8) .GE. 4
C
      ICOUNT = ICALP
      IF (NCALL .EQ. 2) ICOUNT = ICBET
      IF (ITER .GT. 1) GO TO 140
C
C     ----- ITER = 1 -----
C
      DO 120 I = 1,L2
      H1(I) = H0(I)
      H2(I) = H0(I)
  120 CONTINUE
      GO TO 660
C
C     ----- CURRENT FOCK MATRIX IS IN -H0-
C           GET PREVIOUS FOCK MATRICES (H1,H2,H3) -----
C
  140 CALL RDISK(H1,H2,H3,NDAF,L2)
      IF (ITER .GT. 2) GO TO 160
C
C     ----- ITER = 2 -----
C
      IF (DAMPH) GO TO 320
      IF ( .NOT. VSHIFT) GO TO 420
      DAMP = ONE
      GO TO 320
C
C     ----- ITER > 2 -----
C
  160 IF (VSHIFT) GO TO 220
      IF (DAMPRE) GO TO 200
      IF ( .NOT. DAMPH) GO TO 420
      IF (ABS(DE) .GT. EXTTOL) GO TO 180
      IF (DE .GT. ZERO .AND. ICOUNT .GT. IEXTIN) GO TO 180
      GO TO 420
C
  180 DMPTOL = DMPTOL/SHRNKF
      EXTTOL = EXTTOL/SHRNKF
      GO TO 300
C
  200 IF (ABS(DE) .GT. DMPTOL) GO TO 300
      IF (DE .GT. ZERO) GO TO 300
      IF (DAMP .GT. DMPMIN) GO TO 300
      GO TO 420
C
  220 IF ( .NOT. EXTPRE) GO TO 260
      IF (ABS(DE) .GT. EXTTOL) GO TO 240
      IF (DE .GT. ZERO .AND. ICOUNT .GT. IEXTIN) GO TO 240
      GO TO 420
C
  240 EXTTOL = EXTTOL/SHRNKF
      DMPTOL = DMPTOL/SHRNKF
      RRSHFT = RSMIN
      ITERV = 0
      IF (DAMPH) GO TO 300
      GO TO 620
C
  260 IF (ABS(DE) .GT. DMPTOL) GO TO 280
      IF (DE .GT. ZERO) GO TO 280
      IF (DAMPH .AND. DAMP .GT. DMPMIN) GO TO 280
      IF ( .NOT. EXTRAH) GO TO 280
      IF (RRSHFT .GE. VSHTOL) GO TO 280
      IF (ITERV .EQ. 0) GO TO 280
      GO TO 420
C
  280 IF (DAMPH) GO TO 300
      IF (ITERV .LT. 2) GO TO 620
      IF ( .NOT. EXTRAH) GO TO 620
      GO TO 440
C
C     ----- DAVIDSON'S DAMPING -----
C
  300 IF (ICOUNT .LT. IEXTIN) GO TO 320
      IF ( .NOT. VSHIFT .OR. ITERV .GE. 2) GO TO 360
C
  320 DO 340 I = 1,L2
      H0(I) = (H0(I)+DAMP*H1(I))/(ONE+DAMP)
  340 CONTINUE
      GO TO 400
C
C     ----- ITER > 2 , DAMPING -----
C
  360 CUTOFF = PT5-DAMP0
      DO 380 I = 1,L2
      H0(I) = (H0(I)+DAMP*H1(I))/(ONE+DAMP)
  380 CONTINUE
      IF ( .NOT. EXTRAH .OR. CUTOFF .LT. ZERO) GO TO 400
      DAMPRE = .TRUE.
      EXTPRE = .TRUE.
      GO TO 460
C
  400 IF (NCALL .NE. 1) GO TO 660
      DAMPRE = .TRUE.
      EXTPRE = .FALSE.
      GO TO 660
C
  420 IF ( .NOT. EXTRAH) GO TO 620
C
C     ----- POPLE'S EXTRAPOLATION PROCEDURE
C
      IF (NCALL .NE. 1) GO TO 460
      DAMPRE = .FALSE.
      EXTPRE = .TRUE.
  440 DAMP = ZERO
  460 CALL WDISK(H0,H1,H2,NDAF,LL2)
C
C     ----- SKIP TO END IF FIRST CYCLE OR AFTER EXTRAPOLATION -----
C
      DO 480 I = 1,L2
      H3(I) = H2(I)-H3(I)
      H2(I) = H1(I)-H2(I)
      H1(I) = H0(I)-H1(I)
  480 CONTINUE
      IF (ICOUNT .LT. IEXTIN .OR. ITER .LT. 4) GO TO 680
C
C     ----- FIND DISPLACEMENT DP1,DP2,DP3 -----
C
      IF (ITYP .EQ. 2) GO TO 500
      SP11 = TRACEP(H1,H1,L1)
      SP12 = TRACEP(H2,H1,L1)
      SP13 = TRACEP(H3,H1,L1)
      SP22 = TRACEP(H2,H2,L1)
      SP23 = TRACEP(H3,H2,L1)
      SP33 = TRACEP(H3,H3,L1)
      GO TO 520
C
  500 CONTINUE
      SP11 = DDOT(L2,H1,1,H1,1)
      SP12 = DDOT(L2,H2,1,H1,1)
      SP13 = DDOT(L2,H3,1,H1,1)
      SP22 = DDOT(L2,H2,1,H2,1)
      SP23 = DDOT(L2,H3,1,H2,1)
      SP33 = DDOT(L2,H3,1,H3,1)
  520 CONTINUE
      DP1 = SQRT(SP11)
      DP2 = SQRT(SP22)
      DP3 = SQRT(SP33)
C
C        NEXT PROTECT AGAINST SIMPLISTIC EXTRAPOLATION RESTARTING
C        BY MAKING THE "HISTORICAL" FOCK MATRICES ALL IDENTICAL.
C        IN THAT CASE, LOOP 480 HAS MADE ALL DOT PRODUCTS ABOVE ZERO.
C
      IF(DP1.LT.1.0D-10) GO TO 680
      IF(DP2.LT.1.0D-10) GO TO 680
      IF(DP3.LT.1.0D-10) GO TO 680
C
C     ----- FIND COSINE OF ANGLE BETWEEN SUCCESSIVE DISPLACEMENTS -----
C
      COSPHI = SP12/(DP1*DP2)
C
C     ----- FIND COSINE OF ANGLE BETWEEN -DP(3)- AND
C           PLANE OF =DP(1)- AND -DP(2)-.
C
      Z = SP11*SP22-SP12*SP12
      IF ( ABS(Z) .LT. TM17 ) GO TO 680
      X = (SP13*SP22-SP12*SP23)/Z
      Y = (SP23*SP11-SP12*SP13)/Z
      COSPSI = SQRT(X*X*SP11+Y*Y*SP22+TWO*X*Y*SP12)/DP3
C
C     ----- DO NOT EXTRAPOLATE UNLESS -4- CONSECUTIVE POINTS ARE
C           NEARLY COPLANAR -----
C
      IF (COSPSI .LE. TOL) GO TO 680
      IF (DAMP .GT. DMPMIN) GO TO 680
C
C     ----- EXPRESS -DP(1)- AS X*DP(3)(PROJECTED)+Y*DP(2) -----
C
      Y = -Y/X
      X = ONE/X
C
C     ----- TEST IF 2*2 MATRIX HAS REAL EIGENVALUES
C           BETWEEN -TOL/2 AND +TOL/2 -----
C
      XY = Y*Y+FOUR*X
      IF (XY .LT. ZERO) GO TO 680
      XY = ABS(Y)+SQRT(XY)
      IF (XY .LE. TOL1) GO TO 560
C
C     ----- IF -4- POINT EXTRAPOLATION IS NOT POSSIBLE,
C           TRY -3- POINT
C
      IF (ABS(COSPHI) .LE. TOL2) GO TO 680
      X = DP1/(DP2*COSPHI-DP1)
      DO 540 I = 1,L2
      H0(I) = H0(I)+X*H1(I)
  540 CONTINUE
      GO TO 600
C
  560 XXX = X/(ONE-X-Y)
      YYY = (X+Y)/(ONE-X-Y)
      DO 580 I = 1,L2
      H0(I) = H0(I)+XXX*H2(I)+YYY*H1(I)
  580 CONTINUE
C
  600 CONTINUE
      IF(NCALL.EQ.1) ICALP=0
      IF(NCALL.EQ.2) ICBET=0
      CALL DAWRIT(IDAF,IODA,H0,LL2,NDAF,0)
      GO TO 680
C
C     ----- NO DAMPING OR EXTRAPOLATION -----
C
  620 IF (NCALL .NE. 1) GO TO 640
      DAMPRE = .FALSE.
      EXTPRE = .FALSE.
  640 DAMP = ZERO
C
C     ----- SAVE NEW (MODIFIED) FOCK MATRIX -----
C
  660 CALL WDISK(H0,H1,H2,NDAF,LL2)
C
  680 RETURN
      END
C*MODULE SCFLIB  *DECK LEVELS
      SUBROUTINE LEVELS(H,S,V,SHFTMO,SHFTAO,SV,WRK,L0,L1,L2,L3,
     *                  ITER,NDOCC,NHOCC,DE,DEP,DIFF,ITERV,
     *                  NCALL,EXTRAH,CVGING)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION H(L2),S(L2),V(L3),SHFTMO(L2),SHFTAO(L2),SV(L3),WRK(L1)
C
      LOGICAL EXTRAH,CVGING
C
      SAVE DIFFP,PSHIFT
C
      COMMON /ACONV / RRSHFT,EXTTOL,DMPTOL,VSHTOL,IEXTIN
C
      PARAMETER (ZERO=0.0D+00,  HALF=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
C     ----- SHIFT THE DIAGONAL ELEMENTS OF THE MO FOCK MATRIX -----
C     THE SHIFT IS BASED ON CHANGES IN THE ENERGY AND DENSITY
C     THE SHIFT IS CONVERTED TO THE AO BASIS, AND SUMMED INTO -H-
C     THIS ROUTINE REQUIRES -H-, -S-, AND -V- ON ENTRY.
C     -SHFTMO-, -WRK-, -SHFTAO-, AND -SV- ARE ALL SCRATCH STORAGE.
C
      IF(ITER.EQ.1) THEN
         ITERV = 1
         PSHIFT = RRSHFT
         DIFFP = DIFF
         RETURN
      END IF
C
C        NCALL IS USED TO APPLY IDENTICAL SHIFT TO BETA IN UHF
C        NCALL=1 FOR FIRST OPERATOR ONLY
C
      IF(NCALL.NE.1) GO TO 140
C
      ITERV = ITERV+1
      PSHIFT = RRSHFT
      DIFFP = DIFF
C
C     ----- DECREASE THE LEVEL AFTER NSAME ITERATIONS -----
C
      IPWR = (ITER-2)/NSAME
      RSHC = RBASE*RFACT**IPWR
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)
C
C     ----- LET THE LEVEL DROP AFTER NSAME ITERATIONS -----
C
      IF (MOD(ITERV,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 -----
C     ----- CHECK FOR SWITCH TO EXTRAPOLATION ONLY -----
C
      IF(CVGING) RRSHFT = ZERO
      IF(EXTRAH) THEN
         IF(PSHIFT.EQ.ZERO   .AND. ABS(DE).LT.EXTTOL) RRSHFT = ZERO
         IF(RRSHFT.LT.VSHTOL .AND. ABS(DE).LT.DMPTOL) RRSHFT = ZERO
      END IF
      IF(RRSHFT .NE. PSHIFT) ITERV = 0
C
C     ----- GENERATE THE LEVEL SHIFT TO THE FOCK MATRIX -----
C
  140 CONTINUE
      IF(RRSHFT.EQ.ZERO) RETURN
      CALL VCLR(SHFTMO,1,L2)
C
C     ----- SHIFT ANY HALF FILLED ORBITALS BY HALF OF RRSHFT -----
C
      IF(NHOCC.GT.0) THEN
         IS = NDOCC+1
         IE = NDOCC+NHOCC
         DO 160 I = IS,IE
            II = (I*I+I)/2
            SHFTMO(II) = HALF*RRSHFT
  160    CONTINUE
      END IF
C
C     ----- SHIFT ANY VIRTUAL ORBITALS BY RRSHFT -----
C
      IVIRT = NDOCC+NHOCC+1
      DO 200 I = IVIRT,L1
         II = (I*I+I)/2
         SHFTMO(II) = RRSHFT
  200 CONTINUE
C
C     ----- CONVERT THE SHIFT TO THE AO BASIS, AND APPLY IT -----
C
      CALL TFTRIB(SHFTAO,SHFTMO,S,V,SV,WRK,L0,L1,L2,L3)
      CALL VADD(SHFTAO,1,H,1,H,1,L2)
      RETURN
      END
C*MODULE SCFLIB  *DECK MAKEIJ
      SUBROUTINE MAKEIJ(V,FV,E,T,M,M0,N,N0,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION FV(NDIM,M),V(NDIM,M0),E(NDIM,M),T(M0)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- E = V * FV -----
C
      DO 140 J = 1,M
         DO 120 I = 1,M0
            DUM = ZERO
            DO 100 K = 1,N
               DUM = DUM+V(K,I)*FV(K,J)
  100       CONTINUE
            T(I) = DUM
  120    CONTINUE
         DO 130 I = 1,M0
            DUM = T(I)
            IF (J.LE.N0) DUM = DUM+DUM
            E(I,J) = DUM
  130    CONTINUE
  140 CONTINUE
      RETURN
      END
C*MODULE SCFLIB  *DECK MAKFV
      SUBROUTINE MAKFV(F,V,FV,T,IA,M1,M2,N,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION FV(NDIM,*),V(NDIM,*),F(*),T(*),IA(*)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- FV = F * V -----
C
      DO 140 J = M1,M2
         DO 120 I = 1,N
            DUM = ZERO
            DO 100 K = 1,N
               II = MAX(I,K)
               KK = MIN(I,K)
               IK = IA(II)+KK
               DUM = DUM+F(IK)*V(K,J)
  100       CONTINUE
            T(I) = DUM
  120    CONTINUE
         DO 130 I = 1,N
            FV(I,J) = T(I)
  130    CONTINUE
  140 CONTINUE
      RETURN
      END
C*MODULE SCFLIB  *DECK MVOS
      SUBROUTINE MVOS(VEC,EIG,FAO,FMO,SCR,IWRK,WRK,D,BUFFX,BUFFI,
     *                DSH,GHONDO,DDIJ,XINTS,DIRSCF,SCHWRZ,NOPK,NINTMX,
     *                INTTYP,L0,L1,L2,L3,NSH2,MAXG,NFROZ,MVOQ,NE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      LOGICAL DIRSCF,SCHWRZ,GOPARR,DSKWRK,MASWRK,SVDSKW
C
      DIMENSION VEC(L1,L1),EIG(L1),FAO(L2),FMO(L2),SCR(L1,8),
     *          IWRK(L1),WRK(L1,L1),D(L2),BUFFX(NINTMX),BUFFI(NINTMX),
     *          DSH(NSH2),GHONDO(MAXG),DDIJ(*),XINTS(NSH2)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NEX,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (HALF=0.5D+00, ONE=1.0D+00, TWO=2.0D+00, ZERO=0.0D+00)
      PARAMETER (SMALL=1.0D-10)
C
      DATA RHF/8HRHF     /, ROHF/8HROHF    /
C
C     ----- CONSTRUCT MODIFIED VIRTUAL ORBITALS -----
C     C.W.BAUSCHLICHER, JR.  J.CHEM.PHYS.  72,880-885(1980).
C
C     THE FIRST -NFROZ- ORBITALS IN VEC AND EIG WILL BE UNALTERED.
C     THE RHF FOCK MATRIX FOR A CATION OF CHARGE MVOQ (MOVQ+1,
C     IF NEED BE FOR A CLOSED SHELL) WILL BE CONSTRUCTED.  THE
C     CATION FOCK OPERATOR WILL BE TRANSFORMED TO THE BASIS
C     SPANNED BY THE VIRTUAL ORBITALS, DIAGONALIZED, AND THE
C     RESULTING VECTORS AND VALUES USED TO REPLACE THE ORIGINAL
C     VIRTUAL ORBITALS.  TYPICAL CHARGE -MVOQ- SHOULD BE +6 OR SO.
C
      IF(MVOQ.GT.0) THEN
         NOCC = (NE-MVOQ)/2
         NGONE = NE - 2*NOCC
         NVIRT = L0 - NFROZ
         LVIRT = NFROZ+1
         IF (MASWRK) WRITE(IW,9000) NVIRT,NFROZ,NGONE
         IF(SCFTYP.EQ.ROHF  .AND.  NOCC.GT.NB) WRITE(IW,*)
     *     'WARNING: MVOQ DOESN''T REMOVE ALL SINGLY OCCUPIED ELECTRONS'
      ELSE
         NCORE = NUMCOR()
         NOCC  = NFROZ
         NVAL  = NFROZ - NCORE
         NVIRT = L0 - NFROZ
         LVIRT = NFROZ+1
         IF(MASWRK) WRITE(IW,9010) NVIRT,NFROZ,NVAL
      END IF
C
      CALL DAREAD(IDAF,IODA,VEC,L3,15,0)
      CALL DAREAD(IDAF,IODA,EIG,L1,17,0)
C
C    ----- REORDER VIRTUALS BY SYMMETRY GROUPS
C    ----- THIS ALLOWS THE JACOBI DIAGONALIZATION TO MORE EASILY
C    ----- PRESERVER SYMMETRY LATER ON
C
      CALL MVOORD(VEC,EIG,SCR,IWRK,FMO,FAO,WRK,D,L0,L1,L2,L3,NFROZ)
C
      DO 120 I=1,L1
         IWRK(I) = (I*I-I)/2
  120 CONTINUE
C
C     ----- FORM DENSITY MATRIX FOR THE CATION -----
C
      CALL VCLR(SCR,1,L1)
      IF(MVOQ.GT.0) THEN
         DO I=1,NOCC
            SCR(I,1)=TWO
         ENDDO
      ELSE
C            THIS HALF-FILLS RHF OR ROHF VALENCE ORBITALS, BUT NOT GVB
         DO I=1,NCORE
            SCR(I,1)=TWO
         ENDDO
         DO I=NCORE+1,NB
            SCR(I,1)=ONE
         ENDDO
         DO I=NB+1,NA
            SCR(I,1)=ONE/TWO
         ENDDO
      END IF
      CALL DMTX(D,VEC,SCR,NOCC,L1,L1)
C
C     ----- CONSTRUCT SKELETON FOCK MATRIX FOR THE CATION -----
C
      IF(DIRSCF) THEN
         DUMMY = 0.0D+00
         IF(SCHWRZ) CALL SHLDEN(RHF,D,DUMMY,DUMMY,DSH,IWRK,L1,L2,NSH2,1)
         CALL VCLR(FAO,1,L2)
         CALL TWOEI(RHF,DIRSCF,.FALSE.,.FALSE.,.FALSE.,
     *              INTTYP,SCHWRZ,NINT,NSCHWZ,L1,L2,
     *              DUMMY,DUMMY,IDUMMY,1,
     *              XINTS,NSH2,GHONDO,MAXG,DDIJ,
     *              IWRK,D,FAO,DUMMY,DUMMY,DSH,
     *              DUMMY,DUMMY,1)
         CALL DSCAL(L2,HALF,FAO,1)
         II=0
         DO 210 I=1,L1
            II = II+I
            FAO(II) = FAO(II) + FAO(II)
  210    CONTINUE
      ELSE
         SVDSKW = DSKWRK
         DSKWRK = .TRUE.
         CALL SEQREW(IS)
         CALL HSTAR(D,FAO,BUFFX,BUFFI,NINTMX,IWRK,NOPK,.FALSE.)
         DSKWRK = SVDSKW
      END IF
C
C     ----- COMPLETE THE FOCK MATRIX -----
C
      IF (GOPARR) CALL DDI_GSUMF(1004,FAO,L2)
      CALL SYMH(FAO,FMO,IWRK)
      CALL DAREAD(IDAF,IODA,FMO,L2,11,0)
      CALL VADD(FAO,1,FMO,1,FAO,1,L2)
C
C     ----- TRANSFORM TO VIRTUAL MO BASIS, AND MAKE MVO-S -----
C     USE JACOBI DIAGONALIZATION TO TRY TO AVOID MIXING DEGENERATE MO-S
C
      CALL TFTRI(FMO,FAO,VEC(1,LVIRT),SCR,NVIRT,L1,L1)
      DO 300 I=1,(NVIRT*NVIRT+NVIRT)/2
        IF (ABS(FMO(I)) .LT. SMALL) FMO(I)=ZERO
  300 CONTINUE
      IGERR=0
      KDSAVE = KDIAG
      KDIAG = 3
      CALL GLDIAG(L1,NVIRT,NVIRT,FMO,SCR,EIG(LVIRT),
     *            VEC(1,LVIRT),IGERR,IWRK)
      KDIAG = KDSAVE
      IF(IGERR.GT.0) CALL ABRT
      CALL DAREAD(IDAF,IODA,WRK,L3,15,0)
      CALL TFSQB(VEC(1,LVIRT),WRK(1,LVIRT),WRK,NVIRT,L1,L1)
C
      CALL DAWRIT(IDAF,IODA,VEC,L3,15,0)
      CALL DAWRIT(IDAF,IODA,EIG,L1,17,0)
      RETURN
C
 9000 FORMAT(/1X,'GENERATING',I4,' MVO-S (',I4,' ORBITALS FROZEN),',
     *           ' BY REMOVING',I3,' ELECTRONS.')
 9010 FORMAT(/1X,'GENERATING',I4,' MVO-S (',I4,' ORBITALS FROZEN),',
     *           ' BY HALF-FILLING',I4,' VALENCE MOS.')
      END
C
C*MODULE SCFLIB  *DECK MVOORD
      SUBROUTINE MVOORD(VEC,EIG,SCR,IWRK,IWRK2,ISYM,Q,S,
     *                  L0,L1,L2,L3,NFROZ)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION VEC(*),EIG(*),SCR(*),IWRK(*),
     *          ISYM(*),Q(*),S(*),IWRK2(*)
C
      PARAMETER (MXSH=5000, MXATM=2000)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TT(432),INVT(48),NT
C
C     CHECK TO SEE IF WE ARE C1, AND JUST WASTING OUR TIME
C
      IF (NT .EQ. 1) RETURN
C
C     GENERATE SYMMETRY LABELS FOR THE ORBITALS
C     S IS THE OVERLAP MATRIX
C     Q IS THE Q MATRIX
C     ISYM IS THE SYMMETRY LABLES
C
      CALL DAREAD(IDAF,IODA,S,L2,12,0)
      CALL DAREAD(IDAF,IODA,Q,L3,45,0)
      CALL SYMMOS(ISYM,Q,S,VEC,SCR,L0,L1,L0,L1)
      IEND=NFROZ+1
C
C     LOOP OVER EACH IRREP AND ORGANIZE THEM ONE AT A TIME
C
 100  CONTINUE
      ICOUNT=IEND
      IDUMSYM=ISYM(IEND)
      DO 111 I=1,IEND-1
        IWRK(I)=I
 111  CONTINUE
C     MOVE THEM TO THE BEGINING
      DO 115 I=IEND,L0
        IF(ISYM(I) .EQ. IDUMSYM) THEN
          IWRK(ICOUNT)=I
          ICOUNT=ICOUNT+1
        END IF
 115  CONTINUE
C     MOVE OTHER IRREPS TO THE END
      IDUM=ICOUNT
      DO 113 I=IEND+1,L0
        IF(ISYM(I) .NE. IDUMSYM) THEN
          IWRK(IDUM)=I
          IDUM=IDUM+1
        END IF
 113  CONTINUE
C
C     WE NOW HAVE AN ARRAY (IWRK) THAT TELLS HOW TO REORDER THE
C     ORBITALS TO MOVE ALL OF THE ORBITALS OF SYMMETRY TYPE
C     IDUMSYM TO THE BEGINNING.
C
C     WE NOW REORDER VEC, EIG, AND ISYM USING THIS INFORMATION
C
      DO 300 I=1,L0
         IWRK2(I)=IWRK(I)
 300  CONTINUE
      CALL REORDR(VEC,IWRK2,L0,L1)
      DO 310 I=1,L0
         IWRK2(I)=IWRK(I)
 310  CONTINUE
      CALL REORDR(EIG,IWRK2,L0,1)
C     THIS DESTROYS IWRK, CANNOT USE REORDR SINCE ISYM IS INTEGERS
      DO 240 I = 1,L0
         J = IWRK(I)
         ITEMP   = ISYM(I)
         ISYM(I) = ISYM(J)
         ISYM(J) = ITEMP
         DO 220 K = I+1,L0
            IF(IWRK(K).EQ.I) IWRK(K) = J
  220    CONTINUE
  240 CONTINUE
      IEND=ICOUNT
C
C     IF WE ARE NOT FINISHED, GO BACK AND GET THE NEXT SYMMETRY TYPE
C
      IF(IEND .LT. L0) GO TO 100
C
C     WE ARE DONE, SO WRITE THE REORDERED VEC OUT TO DISK
C
      CALL DAWRIT(IDAF,IODA,VEC,L3,15,0)
      RETURN
      END
C*MODULE SCFLIB  *DECK ORBKIN
      SUBROUTINE ORBKIN
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXRT=100)
C
      LOGICAL BETA
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FMCOM / XX(1)
      COMMON /FUNCT / E,EG(3*MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
      DATA UHF,ROHF/8HUHF     ,8HROHF    /
C
      IF(MASWRK) WRITE(IW,9000)
      BETA = SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1 * L1
C
      LOADFM=0
      CALL VALFM(LOADFM)
      I10 = 1 + LOADFM
      I20 = I10 + L3
      I30 = I20 + L1
      I40 = I30 + L2
      LAST = I40 + L2
      NEED = LAST - I10
      CALL GETFM(NEED)
C
      CALL DAREAD(IDAF,IODA,XX(I40),L2,13,0)
      ESUM=ZERO
C
      IF(MASWRK) WRITE(IW,9005)
      CALL DAREAD(IDAF,IODA,XX(I10),L3,15,0)
      DO 20 KK=1,NA
         DO 10 I=1,L1
            XX(I20+I-1)=ZERO
            IF(I.EQ.KK) XX(I20+I-1)=ONE
   10    CONTINUE
         CALL DMTX(XX(I30),XX(I10),XX(I20),NA,L1,L1)
         ETMP=TRACEP(XX(I30),XX(I40),L1)
         IF(MASWRK) WRITE(IW,9010) KK,ETMP
         ESUM=ESUM+ETMP
   20 CONTINUE
C
      IF(BETA) THEN
         IF(MASWRK) WRITE(IW,9006)
         CALL DAREAD(IDAF,IODA,XX(I10),L3,19,0)
         DO 40 KK=1,NB
            DO 30 I=1,L1
               XX(I20+I-1)=ZERO
               IF(I.EQ.KK)XX(I20+I-1)=ONE
   30       CONTINUE
            CALL DMTX(XX(I30),XX(I10),XX(I20),NA,L1,L1)
            ETMP=TRACEP(XX(I30),XX(I40),L1)
            IF(MASWRK) WRITE(IW,9010) KK,ETMP
            ESUM=ESUM+ETMP
   40    CONTINUE
      ELSE
         ESUM=ESUM*TWO
      END IF
C
      IF(MASWRK  .AND.  ABS(ESUM-EKIN).GT.1.0D-06)
     *       WRITE(IW,9020) ESUM,EKIN
      IF(MASWRK) WRITE(IW,9025) ETOT+ESUM
      CALL DAREAD(IDAF,IODA,C,3*NAT,1,0)
      GVIR = DDOT(3*NAT,C,1,EG,1)
      IF(MASWRK) WRITE(IW,9030) ETOT,EKIN,EPOT,GVIR
      CALL RETFM(NEED)
      RETURN
C
 9000 FORMAT(/10X,15(1H-)/10X,15HVIRIAL ANALYSIS/10X,15(1H-))
 9005 FORMAT(/10X,'ALPHA ORBITALS')
 9006 FORMAT(/10X,'BETA  ORBITALS')
 9010 FORMAT(5X,'KINETIC ENERGY OF ORBITAL',I4,' IS',F13.8)
 9020 FORMAT(/,5X,'WARNING! SUM OF ORBITAL KINETIC ENERGIES =',F13.8,
     *       ' DEVIATES SIGNIFICANTLY FROM EKIN =',F13.8)
 9025 FORMAT(/,5X,'SUM OF ETOT AND EKIN',14X,'=',
     *        F19.8)
 9030 FORMAT(/,5X,'TOTAL ENERGY',22X,'=',F19.8,/,
     *       5X,'TOTAL KINETIC ENERGY',14X,'=',F19.8,/,
     *       5X,'ELECTRONIC POTENTIAL ENERGY',7X,'=',F19.8,/,
     *       5X,'POTENTIAL ENERGY FROM RDE/DR TERM =',F19.8,/)
      END
C
C*MODULE SCFLIB  *DECK OVLSEL
      SUBROUTINE OVLSEL(VNEW,ENEW,VOLD,S,SMO,WRK,IWRK,IWRK2,L0,L1,L2)
C
C         NOTE THAT WRK AND IWRK2 MAY BE EQUIVALENT BY CALL
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL MFRZ,RSTRCT
C
      DIMENSION VNEW(L1,L0),ENEW(L0),VOLD(L1,L0),S(L2),
     *          SMO(L1,L0),WRK(L0),IWRK(L0),IWRK2(L0)
C
      PARAMETER (MXAO=8192)
C
      COMMON /MFRPAR/ MFRZ,NUMFRZ,IFRZ(MXAO)
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
      PARAMETER (ZERO=0.0D+00)
C
C        ----- RESTRICT ORBITAL INTERCHANGES DURING SCF -----
C        -VNEW- AND -ENEW- MAY POSSIBLY BE PERMUTED ON EXIT,
C        WHILE -VOLD- AND -S- ARE UNCHANGED ON EXIT.
C
C        ----- CALCULATE OVERLAP BETWEEN THE TWO ORBITAL SETS -----
C        SMO = VOLD-DAGGER * S * VNEW
C
      CALL MTARBR(S,L1,VNEW,L0,SMO,L1,1)
C
      DO 190 J=1,L0
         DO 170 I=1,L0
            DUM = ZERO
            DO 160 K=1,L1
               DUM = DUM + VOLD(K,I)*SMO(K,J)
  160       CONTINUE
            WRK(I) = DUM
  170    CONTINUE
         DO 180 I=1,L0
            SMO(I,J) = WRK(I)
  180    CONTINUE
  190 CONTINUE
C
C        ----- OBTAIN PERMUTATION TO MAXIMUM COINCIDENCE -----
C        THE -GVB- CODE LIKES PHASES TO BE CONSISTENT AS WELL
C
      DO 290 I=1,L0
         SMAX = ZERO
         IWRK(I) = 0
         DO 230 J=1,L0
            DO 210 K=1,I
               IF(IWRK(K).EQ.J) GO TO 230
  210       CONTINUE
            SS = ABS(SMO(I,J))
            IF(SS.GT.SMAX) THEN
               SMAX = SS
               IWRK(I) = J
            END IF
  230    CONTINUE
         IF(SMO(I,IWRK(I)).LT.ZERO) THEN
            DO 260 J=1,L1
               VNEW(J,IWRK(I)) = -VNEW(J,IWRK(I))
  260       CONTINUE
         END IF
  290 CONTINUE
C
C        ----- AND APPLY THAT PERMUTATION -----
C
C     IF MFRZ, ONLY THE FIRST NUMFRZ MO-S ARE RESTRICTED.
C     HOWEVER, IF RESTRICT IS ALSO ON, ALL ORBITALS ARE RESTRICTED.
C
      RSTRCT = MOD(MCONV,16) .GE. 8
C
      IF(MFRZ .AND. .NOT.RSTRCT)THEN
        II = 0
        DO 400 I=NUMFRZ+1,L0
  295      II=II+1
           DO 300 J=1,NUMFRZ
              IF(IWRK(J).EQ.II) GO TO 295
  300      CONTINUE
        IWRK(I)=II
  400   CONTINUE
      END IF
      CALL ICOPY(L1,IWRK,1,IWRK2,1)
      CALL REORDR(VNEW,IWRK,L0,L1)
      CALL ICOPY(L1,IWRK2,1,IWRK,1)
      CALL REORDR(ENEW,IWRK,L0,1)
      RETURN
      END
C*MODULE SCFLIB  *DECK RDISK
      SUBROUTINE RDISK(A1,A2,A3,NDAF,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      DIMENSION A1(2),A2(2),A3(2)
C
      CALL DAREAD(IDAF,IODA,A1,L2,NDAF  ,0)
      CALL DAREAD(IDAF,IODA,A2,L2,NDAF+1,0)
      CALL DAREAD(IDAF,IODA,A3,L2,NDAF+2,0)
C
      RETURN
      END
C*MODULE SCFLIB  *DECK RHR
      SUBROUTINE RHR
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      COMMON /HSYM  / T(35,35),MINK,MAXK,LKT,MINL,MAXL,LLT,NTR
C
      DIMENSION V(35)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- RIGHT MULTIPLY  T  BY  R,
C           RESULT BACK IN  T
C
      GO TO (500,400,300,200,100),LLT
C
C     ----- G SHELL
C
  100 NG=15*(NTR-1)-20
      DO 130 K=MINK,MAXK
      DO 120 L=21,35
      DUM= ZERO
      DO 110 N=21,35
  110 DUM=DUM+T(K,N)*GTR(N-20,NG+L)
  120 V(L)=DUM
      DO 130 L=21,35
  130 T(K,L)=V(L)
      GO TO 500
C
C     ----- F SHELL
C
  200 NF=10*(NTR-1)-10
      DO 230 K=MINK,MAXK
      DO 220 L=11,20
      DUM= ZERO
      DO 210 N=11,20
  210 DUM=DUM+T(K,N)*FTR(N-10,NF+L)
  220 V(L)=DUM
      DO 230 L=11,20
  230 T(K,L)=V(L)
      GO TO 500
C
C     ----- D SHELL
C
  300 ND = 6*(NTR-1)- 4
      DO 330 K = MINK,MAXK
      DO 320 L = 5,10
      DUM = ZERO
      DO 310 N = 5,10
  310 DUM = DUM+T(K,N)*DTR(N-4,ND+L)
  320 V(L) = DUM
      DO 330 L = 5,10
  330 T(K,L) = V(L)
      GO TO 500
C
C     ----- P SHELL
C
  400 NP = 3*(NTR-1)- 1
      DO 430 K = MINK,MAXK
      DO 420 L = 2,4
      DUM = ZERO
      DO 410 N = 2,4
  410 DUM = DUM+T(K,N)*PTR(N-1,NP+L)
  420 V(L) = DUM
      DO 430 L = 2,4
  430 T(K,L) = V(L)
  500 CONTINUE
C
C     ----- LEFT MULTIPLY  T  BY R
C           RESULT BACK IN  T
C
      GO TO (1000,900,800,700,600),LKT
C
C     ----- G SHELL
C
  600 NG=15*(NTR-1)-20
      DO 630 L=MINL,MAXL
      DO 620 K=21,35
      DUM=ZERO
      DO 610 N=21,35
  610 DUM=DUM+GTR(N-20,NG+K)*T(N,L)
  620 V(K)=DUM
      DO 630 K=21,35
  630 T(K,L)=V(K)
      GO TO 1000
C
C     ----- F SHELL
C
  700 NF=10*(NTR-1)-10
      DO 730 L=MINL,MAXL
      DO 720 K=11,20
      DUM= ZERO
      DO 710 N=11,20
  710 DUM=DUM+FTR(N-10,NF+K)*T(N,L)
  720 V(K)=DUM
      DO 730 K=11,20
  730 T(K,L)=V(K)
      GO TO 1000
C
C     ----- D SHELL
C
  800 ND = 6*(NTR-1)-4
      DO 830 L = MINL,MAXL
      DO 820 K = 5,10
      DUM = ZERO
      DO 810 N = 5,10
  810 DUM = DUM+DTR(N-4,ND+K)*T(N,L)
  820 V(K) = DUM
      DO 830 K = 5,10
  830 T(K,L) = V(K)
      GO TO 1000
C
C     ----- P SHELL
C
  900 NP = 3*(NTR-1)- 1
      DO 930 L = MINL,MAXL
      DO 920 K = 2,4
      DUM = ZERO
      DO 910 N = 2,4
  910 DUM = DUM+PTR(N-1,NP+K)*T(N,L)
  920 V(K) = DUM
      DO 930 K = 2,4
  930 T(K,L) = V(K)
 1000 CONTINUE
      RETURN
      END
C*MODULE SCFLIB  *DECK SCFIN
      SUBROUTINE SCFIN
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DIIS,DAMPH,EXTRAP,RSTRCT,VSHIFT,COUPLE,UHFNOS,DEM,
     *        DIRSCF,FDIFF,GOPARR,DSKWRK,MASWRK,VTSCAL,VIROK,
     *        LVCLN,SOSCF,NOCONV,ABEL,ABELPT,ACAVO,VVOS,SG1
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXAO=8192)
      PARAMETER (MXATM=2000, MXGRID=10)
      PARAMETER (ZERO=0.0D+00, HALF=0.5D+00, ONE=1.0D+00, TWO=2.0D+00)
C
      DIMENSION KCORB(2,12)
C
C         THE COMMON BLOCKS RELATED TO SCF CONVERGENCE THAT NEED
C         TO BE CONSOLIDATED TO JUST ONE NAMED COMMON BLOCK ARE:
C         ACONV, CONV, DEMOPT, DIISSO, DMPING, OPTSCF, SCFOPT, SCFWFN!
C
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
      COMMON /DEMOPT/ DEMCUT,IDEM
      COMMON /DMPING/ SHIFTO,SHIFTV,DMPCUT,SWDIIS,DIRTHR
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      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 /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /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),MIN(MXSH),MAX(MXSH),NSHELL
      COMMON /OPTSCF/ DIRSCF,FDIFF
      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 /SCFWFN/ AROHF(3),BROHF(3),PACAVO(6),IACAVO,IUHFNO,
     *                MVOQ
      COMMON /VIRIAL/ SCALTE,SCALTT,GVIR,VTCONV,MAXVT,VTSCAL,VIROK,LVCLN
      COMMON /VVOPAR/ IVVOS,NVVOS,BNDDEN
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
C     ----- SET UP NAMELIST $SCF TABLES -----
C     -NCONV- BECAME AN UNDOCUMENTED OPTION IN JUNE 2001, IT WAS
C     REPLACED BY THE FLOATING POINT VALUE -CONV- CARRIED AS -CONVHF-.
C     IN A YEAR OR TWO, THIS FIRST ELEMENT CAN BE RECYLED IF A NEW
C     INPUT VARIABLE IS NEEDED IN THIS GROUP.
C
      PARAMETER (NNAM=39)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA SCF/8HSCF     /
      DATA QNAM/8HNCONV   ,8HNPUNCH  ,8HNCO     ,
     *          8HNSETO   ,8HDMPCUT  ,8HNPAIR   ,8HNO      ,
     *          8HCICOEF  ,8HF       ,8HALPHA   ,8HBETA    ,
     *          8HDIIS    ,8HETHRSH  ,
     *          8HMAXDII  ,8HEXTRAP  ,8HDAMP    ,
     *          8HSHIFT   ,8HRSTRCT  ,8HCOUPLE  ,8HUHFNOS  ,
     *          8HMVOQ    ,8HDEM     ,8HDEMCUT  ,8HDIRSCF  ,
     *          8HFDIFF   ,8HVTSCAL  ,8HVTCONV  ,8HMAXVT   ,
     *          8HSCALF   ,8HSOSCF   ,8HSOGTOL  ,8HCONV    ,
     *          8HACAVO   ,8HPACAVO  ,8HVVOS    ,8HNOCONV  ,
     *          8HDIRTHR  ,8HSWDIIS  ,8HNPREO   /
      DATA KQNAM/1,1,1,   1,3,1,101,   253,253,3253,3253,
     *           0,3,   1,0,0,   0,0,0,0,   1,0,3,0,
     *           0,0,3,1,   3,0,3,3,   0,63,0,0,   3,3,41/
C
      DATA RHF,UHF /8HRHF     ,8HUHF     /
      DATA GVB,ROHF/8HGVB     ,8HROHF    /
      DATA RNONE/8HNONE     /
      DATA NONE/4HNONE/
      DATA BNDANA/8HBONDANAL/
C
C     ----- READ IN THE INPUT DATA NECESARY FOR THE SCF MODULES -----
C
      L1 = NUM
      ABEL = ABELPT()
      CALL DERCHK(NDER)
C
C     ----- SET UP SOME UPPER LIMITS -----
C
      NHAMX = 25
      NSETMX = 10
      NPAIRX = 12
C
C     ----- SET UP TRIANGULAR INDEX MATRIX -----
C
      DO 100 I = 1,L1
         IA(I) = (I*I-I)/2
  100 CONTINUE
C
C        IF THE USER WANTS DIIS AND SOSCF TO BOTH BE OFF, THEY
C        MUST SELECT NOCONV=.TRUE.  OTHERWISE ONE OF THESE WILL
C        BE SWITCHED ON AFTER THE NAMEIO CALL, SEE BELOW.
C
      SOSCF  = .FALSE.
      DIIS   = .FALSE.
      NOCONV = .FALSE.
C
C        AUXILIARY CONVERGERS, USED MAINLY BEFORE SOSCF/DIIS ENGAGE
C
      EXTRAP = .TRUE.
      DAMPH  = .FALSE.
      VSHIFT = .FALSE.
      RSTRCT = .FALSE.
      DEM    = .FALSE.
C                  NUMERICAL VALUES ASSOCIATED WITH CONVERGERS
      DMPCUT = 0.0D+00
      ETHRSH = 0.5D+00
      DEMCUT = 0.5D+00
      SOGTOL = 0.25D+00
      MAXDII = 10
C
C                  POSSIBLE TOLERANCE TO SWITCH FROM DIIS TO SOSCF
C
      SWDIIS = 0.0D+00
C
C     FORCE RECOMPUTATION OF THE FOCK MATRIX(CES) FOR FDIFF IF THE
C     CONVERGENCE STAYS BELOW DIRTHR FOR 5 ITERATIONS.
C
      DIRTHR=ZERO
C
C        IF NCONV IS IN THE INPUT, IT WILL STILL BE USED
C
      NCONV  = 0
      CONVHF = 0.0D+00
      NPUNCH = 2
      UHFNOS = .FALSE.
      ACAVO  = .FALSE.
      VVOS   = .FALSE.
      MVOQ   = 0
      CALL VCLR(PACAVO,1,6)
      PACAVO(6) = -1.0D+00
C
      DIRSCF = .FALSE.
      FDIFF  = .TRUE.
C
      VTSCAL = .FALSE.
      VTCONV = 1.0D-06
      MAXVT  = 20
      SCALF  = 1.0D+00
      VIROK  = .FALSE.
C                      NEXT ONES APPLY ONLY TO -GVB- CASE
      NCO   = 0
      NSETO = 0
      NPAIR = 0
      DO 120 I = 1,NSETMX
         NO(I) = 0
  120 CONTINUE
      DO 130 I=1,NPAIRX
         CICOEF(1,I) = 0.90D+00
         CICOEF(2,I) =-0.20D+00
  130 CONTINUE
      COUPLE=.FALSE.
      N2 = (NHAMX*NHAMX+NHAMX)/2
      CALL VCLR(F,1,NHAMX)
      CALL VCLR(ALPHA,1,N2)
      CALL VCLR(BETA ,1,N2)
      NPREO(1)=1
      NPREO(2)=L1
      NPREO(3)=2
      NPREO(4)=1
C
C     ----- READ IN $SCF NAMELIST -----
C
      JRET = 0
      CALL NAMEIO(IR,JRET,SCF,NNAM,QNAM,KQNAM,
     *            NCONV,NPUNCH,NCO,NSETO,DMPCUT,NPAIR,
     *            NO,CICOEF,F,ALPHA,BETA,DIIS,ETHRSH,MAXDII,
     *            EXTRAP,DAMPH,VSHIFT,RSTRCT,COUPLE,UHFNOS,MVOQ,
     *            DEM,DEMCUT,DIRSCF,FDIFF,VTSCAL,VTCONV,MAXVT,
     *            SCALF,SOSCF,SOGTOL,CONVHF,ACAVO,PACAVO,VVOS,
     *            NOCONV,DIRTHR,SWDIIS,NPREO,      0,0,0,0,0,
     *            0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0)
      IF (JRET .EQ. 2) THEN
         IF (MASWRK) WRITE (IW,9008)
         CALL ABRT
      END IF
      NERR=0
C
C         CHECK VARIOUS ORBITAL MODIFICATION SCHEMES
C
      IF(RUNTYP.EQ.BNDANA) THEN
         VVOS=.TRUE.
         MVOQ=0
         ACAVO=.FALSE.
      END IF
C
      CALL DERCHK(NDER)
      IF((MVOQ.NE.0 .OR. ACAVO  .OR.  VVOS)  .AND.  NDER.GT.0) THEN
         IF(MASWRK) WRITE(IW,*)
     *      '-MVOQ- OR -ACAVO- OR -VVOS- IS INCOMPATIBLE WITH GRADIENTS'
         NERR=NERR+1
      END IF
C
      IF((MVOQ.NE.0 .OR. ACAVO)  .AND.  SCFTYP.EQ.UHF) THEN
         IF(MASWRK) WRITE(IW,*)
     *      'NEITHER -MVOQ- NOR -ACAVO- IS PROGRAMMED FOR SCFTYP=UHF.'
         NERR = NERR+1
      END IF
      IF(MVOQ.LT.0  .AND.  SCFTYP.EQ.GVB) THEN
         IF(MASWRK) WRITE(IW,*)
     *      'MVOQ=-1 IS NOT PROGRAMMED FOR SCFTYP=GVB.'
         NERR = NERR+1
      END IF
      IF(VVOS  .AND.  SCFTYP.NE.RHF) THEN
         IF(MASWRK) WRITE(IW,*)
     *      '-VVOS- IS PROGRAMMED ONLY FOR SCFTYP=RHF.'
         NERR = NERR+1
      END IF
      IF(ACAVO  .AND.  DIRSCF) THEN
         IF(MASWRK) WRITE(IW,*)
     *      '-ACAVO- IS NOT YET PROGRAMMED FOR DIRSCF=.TRUE.'
         NERR = NERR+1
      END IF
C
      NMODVO = 0
      IF(MVOQ.NE.0) NMODVO=NMODVO+1
      IF(ACAVO)     NMODVO=NMODVO+1
      IF(VVOS)      NMODVO=NMODVO+1
      IF(NMODVO.GT.1) THEN
         IF(MASWRK) WRITE(IW,*)
     *      '-MVOQ-, -ACAVO-, AND -VVOS- ARE MUTUALLY EXCLUSIVE'
         NERR = NERR+1
      END IF
      IF(NMODVO.EQ.1) THEN
         IF(MPLEVL.GT.0  .OR.  CCTYP.NE.RNONE) THEN
            IF(MASWRK) WRITE(IW,*)
     *           'MP2 OR CC RUNS REQUIRE CANONICAL VIRTUALS'
            NERR=NERR+1
         END IF
      END IF
      IF(NMODVO.GT.0  .AND.  MPCTYP.NE.NONE) THEN
         IF(MASWRK) WRITE(IW,*)
     *        'MODIFICATION OF VIRTUALS FOR MOPAC IS NOT POSSIBLE'
         NERR=NERR+1
      END IF
C
                 IACAVO=0
      IF(ACAVO)  IACAVO=1
                 IVVOS=0
      IF(VVOS)   IVVOS=1
                 NVVOS=0
                 IUHFNO=0
      IF(UHFNOS) IUHFNO=1
C
C        MOSTLY, SOSCF WILL BE THE DEFAULT CONVERGER, BUT IN
C        A FEW CASES WE EXPECT DIIS WILL BE MUCH BETTER THAN SOSCF
C        OF COURSE, THESE TWO ARE MUTUALLY INCOMPATIBLE.
C
      IF(DIIS .AND. SOSCF) THEN
         NERR = NERR+1
         IF(MASWRK) WRITE(IW,*)
     *      '-SOSCF- AND -DIIS- ARE MUTUALLY EXCLUSIVE CONVERGERS'
      END IF
C
      IF(MPCTYP.NE.NONE) DEM = .FALSE.
      IF(MPCTYP.NE.NONE .AND. .NOT.DIIS .AND. .NOT.SOSCF) NOCONV=.TRUE.
C
C         SOSCF IS USED FOR MOST THINGS, EXCEPT
C               UHF WHICH BASICALLY NEVER CONVERGES WITH SOSCF,
C               DFT WHICH IS A LOT TRICKIER THAN ORDINARY SCF,
C               NON-ABELIAN GROUPS.
C         WHERE EXPERIENCE TEACHES US DIIS IS SINE QUA NON.
C
C         UP TO MARCH 2007, NON-ABELIAN GROUPS IN -ROHF- ONLY
C         WERE SWITCHED TO DIIS, BUT AFTER THAT ALL SCFTYP ARE.
C         DIIS USES SYMDIA TO CONTROL ORBITAL DEGENERACIES,
C         WHEREAS IN SOME RUNS, SOSCF LOSES A BIT OF ORBITAL
C         DEGENERACY AND THUS GETS SLIGHTLY INCORRECT ENERGIES
C
      IF(.NOT.SOSCF  .AND.  .NOT.DIIS  .AND.  .NOT.NOCONV) THEN
         IF(SCFTYP.EQ.UHF  .OR.  .NOT.ABEL
     *                     .OR.  DFTTYP(1).NE.0.0D+00
     *                     .OR.  NDFTFG.NE.0) THEN
            DIIS = .TRUE.
            SOSCF= .FALSE.
         ELSE
            SOSCF  = .TRUE.
            DIIS   = .FALSE.
         END IF
      END IF
C
C         TRIVIAL PROBLEMS SHOULD NOT INVOLVE POWERFUL CONVERGERS
C            FOR FMO NEITHER NA+NB NOR L1 IS MEANINGFUL
      IF(NFG.EQ.0) THEN
         IF(NA+NB.EQ.1) NOCONV=.TRUE.
         IF(L1.LT.10)   NOCONV=.TRUE.
      END IF
C
      IF(NOCONV) DIIS  = .FALSE.
      IF(NOCONV) SOSCF = .FALSE.
C
      MCDEF=0
      IF(EXTRAP) MCDEF=MCDEF+1
      IF(DAMPH)  MCDEF=MCDEF+2
      IF(VSHIFT) MCDEF=MCDEF+4
      IF(RSTRCT) MCDEF=MCDEF+8
      IF(DIIS)   MCDEF=MCDEF+16
      IF(DEM)    MCDEF=MCDEF+32
      IF(SOSCF)  MCDEF=MCDEF+64
      MCONV = MCDEF
C
      IF(.NOT.DIIS.AND.NFG.EQ.0) SWDIIS=0.0D+00
C
C        NEXT CODE ALLOWS OLD INPUT FILES WITH -NCONV- TO WORK.
C        NOTE THAT CORRELATED GRADIENTS REQUIRE A BIT TIGHTER SCF,
C        TO ENSURE VIRTUALS ARE ALSO WELL CONVERGED.
C        DO NOT ALLOW ANYONE TO INPUT A VALUE THAT IS TOO SLOPPY.
C
      IF(CONVHF.NE.ZERO  .AND.  NCONV.NE.0) THEN
         IF(MASWRK) WRITE(IW,9401)
         NERR = NERR+1
      END IF
      IF(CONVHF.EQ.ZERO) THEN
         CONVHF = 1.0D-05
         IF(MPLEVL.GT.0)    CONVHF=1.0D-06
         IF(CCTYP.NE.RNONE) CONVHF=1.0D-06
         IF(CITYP.NE.RNONE) CONVHF=1.0D-06
         IF(NFG.NE.0)       CONVHF=1.0D-07
      END IF
      IF(NCONV.NE.0) CONVHF=10.0D+00**(-NCONV)
      IF(CONVHF.GT.1.0D-03) CONVHF=1.0D-03
C
      IF(MASWRK) THEN
         IF(EXTRAP) WRITE(IW,9402)
         IF(DAMPH)  WRITE(IW,9404)
         IF(VSHIFT) WRITE(IW,9406)
         IF(RSTRCT) WRITE(IW,9408)
         IF(DIIS)   WRITE(IW,9410)
         IF(SOSCF)  WRITE(IW,9412)
         IF(VTSCAL) WRITE(IW,9415) MAXVT,VTCONV
      END IF
C
      IF(DIRSCF  .AND.  DEM) THEN
         IF(MASWRK) WRITE(IW,*) 'DEM AND DIRSCF ARE INCOMPATIBLE.'
         CALL ABRT
      END IF
C
      IF(MASWRK) WRITE(IW,9050) NPREO
C
C     ALLOW INITIAL ORBITAL EXPONENT SCALING FOR VIRIAL RUN
C
      SCALTT = SCALF
      IF (SCALF.NE.1.0D+00) THEN
         IF (MASWRK) WRITE(IW,9416) SCALF
         CALL NORMAO(1)
         DO 90 I=1,MXGTOT
            EX(I)=EX(I)*SCALF
   90    CONTINUE
         CALL NORMAO(2)
      END IF
C
      IF(SCFTYP.EQ.RHF) THEN
         NCO=NA
         NORB=NA
         NSETO=0
         NPAIR=0
         NHAM=1
         ALPHA(1)=TWO
         BETA(1)=-ONE
         F(1)=ONE
      END IF
C
      IF(SCFTYP.NE.GVB  .AND.  SCFTYP.NE.ROHF) THEN
         IF(NERR.GT.0) CALL ABRT
         RETURN
      END IF
C
      IF(SCFTYP.EQ.ROHF) THEN
         NCO   = NB
         NSETO = 0
         IF(NA.GT.NB) NSETO = 1
         NO(1) = NA - NB
         NPAIR = 0
      END IF
C
C     ----- FILL /GVBWFN/ -----
C     COUNT THE NUMBER OF OCCUPIED AND OPEN SHELL ORBITALS
C
      NOPEN = 0
      IF(NSETO.GT.0) THEN
            DO 220 I = 1,NSETO
               NOP = NO(I)
               IF(NOP.LE.0) THEN
                  IF(MASWRK) WRITE(IW,9218) I,NOP
                  NERR = NERR+1
               END IF
               NOPEN = NOPEN+NOP
  220       CONTINUE
         END IF
      NORB = NCO+NOPEN+NPAIR+NPAIR
C
C     ----- CALCULATE THE NUMBER OF FOCK OPERATORS = NHAM -----
C
      NCORES = 0
      IF(NCO.NE.0) NCORES = 1
      NHAM = NCORES+NPAIR+NPAIR+NSETO
C
      IF(SCFTYP.EQ.GVB) THEN
         IF (MASWRK) THEN
         WRITE (IW,9028)
         WRITE (IW,9068) NORB,NCO,NPAIR,NSETO,(NO(I),I = 1,NSETO)
         END IF
         NBASE = NCO
         IF (NSETO.GT.0) THEN
               IF (MASWRK) WRITE (IW,9128)
               DO 300 I = 1,NSETO
                  NOPL = NBASE+1
                  NOPH = NBASE+NO(I)
                  IF (MASWRK) WRITE (IW,9148) I,(J,J = NOPL,NOPH)
                  NBASE = NBASE+NO(I)
  300          CONTINUE
            END IF
         IF (NPAIR.GT.0) THEN
               IF (MASWRK) WRITE (IW,9168)
               DO 340 I = 1,NPAIR
                  NORB1 = NBASE+1
                  NORB2 = NBASE+2
                  IF (MASWRK) WRITE (IW,9188) I,NORB1,NORB2
                  NBASE = NBASE+2
  340          CONTINUE
            END IF
      END IF
C
C     ----- CHECK LIMITS -----
C
      IF (NORB.GT.L1   .OR.  NHAM.GT.NHAMX   .OR.
     *    NSETO.GT.NSETMX   .OR.  NPAIR.GT.NPAIRX) THEN
            IF (MASWRK) WRITE (IW,9208) NORB,L1,NHAM,NHAMX,NSETO,NSETMX,
     *                      NPAIR,NPAIRX
            NERR = NERR+1
         END IF
C
C     ----- GENERATE NCONF -----
C                                     ZERO ALL...
      DO 345 I=1,L1
         NCONF(I) = 0
  345 CONTINUE
C                                   ...THEN SET CORE...
      ICONF = 0
      IF(NCO.GT.0) THEN
         ICONF = ICONF + 1
         DO 350 I=1,NCO
            NCONF(I) = ICONF
  350    CONTINUE
      END IF
C                                   ...THEN ANY OPEN SHELLS...
      LOC = NCO
      IF(NSETO.GT.0) THEN
         DO 360 I=1,NSETO
            NOP = NO(I)
            ICONF = ICONF+1
            DO 355 J=1,NOP
               NCONF(LOC+J) = ICONF
  355       CONTINUE
            LOC = LOC + NOP
  360    CONTINUE
      END IF
C                                    ...AND THEN ANY PAIRS.
      IF(NPAIR.GT.0) THEN
         NP2 = NPAIR+NPAIR
         DO 370 I=1,NP2
            ICONF = ICONF + 1
            NCONF(LOC+I) = ICONF
  370    CONTINUE
      END IF
C
C     ----- SET UP DEFAULT COUPLING CONSTANTS -----
C     NOTE THAT ALL -GVB- CASES ALLOW PAIRS TO BE COUPLED BEHIND
C     THE CLOSED AND OPEN SHELL SPACES, IF DESIRED.
C
      IF(COUPLE) GO TO 500
C
C     ----- HIGH SPIN CASE (ROHF MODULE) -----
C     THE DEFAULT IS THE ROOTHAAN SINGLE MATRIX CASE.
C     FOR OTHER INTERESTING CASES, SEE SUBROUTINE ROFOCK.
C
      IF(SCFTYP.EQ.ROHF) THEN
         ALPHA(1)= -HALF
         ALPHA(2)=  HALF
         ALPHA(3)=  ONE+HALF
         BETA(1) =  ONE+HALF
         BETA(2) =  HALF
         BETA(3) = -HALF
         COUPLE = .TRUE.
      END IF
C
C     ----- CLOSED SHELL SINGLET (GVB MODULE) -----
C
      IF(NSETO.EQ.0  .AND.  MUL.EQ.1  .AND.  SCFTYP.EQ.GVB) THEN
         F(1) = ONE
         ALPHA(1) = TWO
         BETA(1) = -ONE
         COUPLE=.TRUE.
      END IF
C
C     ----- OPEN SHELL DOUBLET (GVB MODULE) -----
C
      IF(NSETO.EQ.1  .AND.  NO(1).EQ.1  .AND.
     *     MUL.EQ.2  .AND.  SCFTYP.EQ.GVB) THEN
         IF(NCO.GT.0) THEN
            F(1) = ONE
            F(2) = HALF
            ALPHA(1) = TWO
            ALPHA(2) = ONE
            ALPHA(3) = HALF
            BETA(1) = -ONE
            BETA(2) = -HALF
            BETA(3) = -HALF
         ELSE
            F(1) = HALF
            ALPHA(1) = HALF
            BETA(1) = -HALF
         END IF
         COUPLE=.TRUE.
      END IF
C
C     ----- OPEN SHELL SINGLET (GVB MODULE) -----
C
      IF(NSETO.EQ.2  .AND.  NO(1).EQ.1  .AND.  NO(2).EQ.1   .AND.
     *     MUL.EQ.1  .AND.  SCFTYP.EQ.GVB) THEN
         IF(NCO.GT.0) THEN
            F(1) = ONE
            F(2) = HALF
            F(3) = HALF
            ALPHA(1) = TWO
            ALPHA(2) = ONE
            ALPHA(3) = ZERO
            ALPHA(4) = ONE
            ALPHA(5) = HALF
            ALPHA(6) = ZERO
            BETA(1) = -ONE
            BETA(2) = -HALF
            BETA(3) =  ZERO
            BETA(4) = -HALF
            BETA(5) =  HALF
            BETA(6) =  ZERO
         ELSE
            F(1) = HALF
            F(2) = HALF
            ALPHA(1) = ZERO
            ALPHA(2) = HALF
            ALPHA(3) = ZERO
            BETA(1) =  ZERO
            BETA(2) =  HALF
            BETA(3) =  ZERO
         END IF
         COUPLE=.TRUE.
      END IF
C
C     ----- OPEN SHELL TRIPLET (GVB MODULE) -----
C
      IF(NSETO.EQ.2  .AND.  NO(1).EQ.1  .AND.  NO(2).EQ.1   .AND.
     *     MUL.EQ.3  .AND.  SCFTYP.EQ.GVB) THEN
         IF(NCO.GT.0) THEN
            F(1) = ONE
            F(2) = HALF
            F(3) = HALF
            ALPHA(1) = TWO
            ALPHA(2) = ONE
            ALPHA(3) = HALF
            ALPHA(4) = ONE
            ALPHA(5) = HALF
            ALPHA(6) = HALF
            BETA(1) = -ONE
            BETA(2) = -HALF
            BETA(3) = -HALF
            BETA(4) = -HALF
            BETA(5) = -HALF
            BETA(6) = -HALF
         ELSE
            F(1) = HALF
            F(2) = HALF
            ALPHA(1) = HALF
            ALPHA(2) = HALF
            ALPHA(3) = HALF
            BETA(1) = -HALF
            BETA(2) = -HALF
            BETA(3) = -HALF
         END IF
         COUPLE=.TRUE.
      END IF
C
C     IF VALUES WERE NOT DEFAULTS AND WERE NOT INPUT, BOMB THE JOB.
C
  500 CONTINUE
      IF(.NOT.COUPLE) THEN
         IF (MASWRK) WRITE (IW,9228) SCFTYP,MUL,NSETO
         NERR = NERR+1
      END IF
C
C     ----- THE ROHF-IN-UHF CODE WORKS WITH AROHF AND BROHF -----
C     WE SET F,ALPHA,BETA FOR THE OPEN SHELL CPHF CODE'S USE.
C
      IF(SCFTYP.EQ.ROHF) THEN
         AROHF(1) = ALPHA(1)
         AROHF(2) = ALPHA(2)
         AROHF(3) = ALPHA(3)
         BROHF(1) =  BETA(1)
         BROHF(2) =  BETA(2)
         BROHF(3) =  BETA(3)
C
         IF(NCO.GT.0) THEN
            IF(NA.EQ.NB) THEN
               F(1) = ONE
               ALPHA(1) = TWO
               BETA(1)  = -ONE
            ELSE
               F(1) = ONE
               F(2) = HALF
               ALPHA(1) = TWO
               ALPHA(2) = ONE
               ALPHA(3) = HALF
               BETA(1) = -ONE
               BETA(2) = -HALF
               BETA(3) = -HALF
            END IF
         ELSE
            F(1) = HALF
            ALPHA(1) = HALF
            BETA(1) = -HALF
         END IF
      END IF
C
      IF(SCFTYP.EQ.ROHF) THEN
         IF(NERR.GT.0) CALL ABRT
         RETURN
      END IF
C
C     ----- MAKE SURE THE $SCF INPUT MATCHES DESIRED CHARGE -----
C
      NELEC = 2*NCO + 2*NPAIR
      ELEC = ZERO
      DO 710 I=1,NSETO
         ELEC = ELEC + TWO * NO(I) * F(I+NCORES)
  710 CONTINUE
      NELEC = NELEC + INT(ELEC+0.0001D+00)
      ICHARG = NE - NELEC + ICH
      IF(ICHARG.NE.ICH) THEN
         IF (MASWRK) WRITE(IW,9070) ICH,ICHARG,NELEC
         NERR = NERR+1
      END IF
      IF(NERR.GT.0) CALL ABRT
C
C     ----- GENERATE KCORB -----
C
      IF(NPAIR.GT.0) THEN
            NBASE = NCO + NOPEN
            DO 760 KPAIR = 1,NPAIR
               KCORB(1,KPAIR) = NBASE + 1
               KCORB(2,KPAIR) = NBASE + 2
               NBASE = NBASE + 2
  760       CONTINUE
         END IF
C
C     ----- NORMALIZE CI COEFFICIENTS, SET PAIR'S F, ALPHA, BETA -----
C
      IF (NPAIR.GT.0) CALL CIEXPR(KCORB,IA,L1)
C
C     ----- PRINT THE COUPLING PARAMETERS -----
C
      IF (MASWRK) THEN
      WRITE (IW,9268)
      WRITE (IW,9288)
      DO 800 I = 1,NHAM
         WRITE (IW,9308) I,F(I)
  800 CONTINUE
      WRITE (IW,9328)
      CALL PRTRI(ALPHA,NHAM)
      WRITE (IW,9348)
      CALL PRTRI(BETA,NHAM)
      IF(NPAIR.GT.0) THEN
            WRITE (IW,9368)
            NBASE = NCO+NOPEN
            DO 820 KPAIR = 1,NPAIR
               WRITE (IW,9388) KPAIR,CICOEF(1,KPAIR),CICOEF(2,KPAIR)
  820       CONTINUE
         END IF
      END IF
      RETURN
C
 9008 FORMAT(1X,'ERROR IN $SCF INPUT - STOP')
 9028 FORMAT(/10X,25(1H*)/10X,'ROHF-GVB INPUT PARAMETERS'/10X,25(1H*)/)
 9050 FORMAT(1X,'ORBITAL PRINTING OPTION: NPREO=',4I4)
 9068 FORMAT(10X,'NORB   =',I5,10X,'NCO    =',I5/
     *       10X,'NPAIR  =',I5,10X,'NSETO  =',I5,:,/
     *       10X,'NO     =',20I5)
 9070 FORMAT(1X,'**** ERROR **** CHECK CHARGE AND $SCF'/
     *       1X,'$CONTRL SAYS ICHARG=',I4,' BUT'/
     *       1X,'$SCF SAYS ICHARG=',I4,' DUE TO',I5,' ELECTRONS.')
 9128 FORMAT(10X,19HOPEN SHELL ORBITALS)
 9148 FORMAT(10X,'SET ',I5,' HAS ORBS',10I5)
 9168 FORMAT(10X,'PAIR ORBITALS')
 9188 FORMAT(10X,'PAIR',I5,' HAS ORBS',2I5)
 9208 FORMAT(1X,'***** ERROR **** DIMENSIONS EXCEEDED *****'/
     *       1X,' NORB=',I5,' NO.AOS=',I5/
     *       1X,' NHAM=',I5,' MAX=',I5/
     *       1X,'NSETO=',I5,' MAX=',I5/
     *       1X,'NPAIR=',I5,' MAX=',I5)
 9218 FORMAT(1X,'***** ERROR *****'/1X,'OPEN SHELL NUMBER',I3,
     *          ' HAS NONSENSE DEGENERACY=',I5)
 9228 FORMAT(//1X,'**** ERROR ****'/1X,'SCFTYP=',A8,' MULT=',I4,
     *            ' NSETO=',I4,' HAS NO DEFAULT COUPLING PARAMETERS'/
     *      1X,'READ THE ROHF-GVB SECTION IN THE REFS.DOC CHAPTER!')
 9268 FORMAT(/10X,28(1H-)/10X,'ROHF-GVB COUPLING PARAMETERS'/
     *        10X,28(1H-)/)
 9288 FORMAT(10X,'F VECTOR (OCCUPANCIES)')
 9308 FORMAT(I5,F15.10)
 9328 FORMAT(10X,'ALPHA COUPLING COEFFICEINTS')
 9348 FORMAT(10X,' BETA COUPLING COEFFICIENTS')
 9368 FORMAT(/10X,28HNATURAL ORBITAL COEFFICIENTS/
     *        10X,'N.O. ',15X,'PAIR CICOEF-S')
 9388 FORMAT(10X,I3,5X,F15.10,2X,F15.10)
 9401 FORMAT(1X,'IT IS ILLOGICAL TO INPUT BOTH NCONV AND CONV IN $SCF')
 9402 FORMAT(10X,'EXTRAPOLATION IN EFFECT')
 9404 FORMAT(10X,'DAMPING IN EFFECT')
 9406 FORMAT(10X,'LEVEL SHIFTING IN EFFECT')
 9408 FORMAT(10X,'RESTRICTION OF ORBITAL MIXING IN EFFECT')
 9410 FORMAT(10X,'DIIS IN EFFECT')
 9412 FORMAT(10X,'SOSCF IN EFFECT')
 9415 FORMAT(10X,'VIRIAL SCALING IN EFFECT, MAXVT=',I4,' VTCONV=',F13.9)
 9416 FORMAT(10X,'INITIAL ORBITAL EXPONENTS WILL BE SCALED BY ',F13.9)
      END
C*MODULE SCFLIB  *DECK SCFSAV
      SUBROUTINE SCFSAV(V,D,E,NDAF,L1,L2,L3)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      DIMENSION D(L2),V(L3),E(L1)
C
      CALL DAWRIT(IDAF,IODA,V,L3,NDAF  ,0)
      CALL DAWRIT(IDAF,IODA,D,L2,NDAF+1,0)
      CALL DAWRIT(IDAF,IODA,E,L1,NDAF+2,0)
C
      RETURN
      END
C*MODULE SCFLIB  *DECK SHLDEN
      SUBROUTINE SHLDEN(SCFTYP,DA,DB,DGVB,DSH,IA,L1,L2,NSH2,NFO)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL ALPHA,BETA,FMOESP
C
      DIMENSION DA(L2),DB(L2),DGVB(L2,NFO),DSH(NSH2),IA(L1)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      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 /FMORUN/ ESPSCF,E0SCF(2),IDAFMO,ICURLAY,ICURFG,JCURFG,
     *                KCURFG,NCURSH,NGAU,ICURPOP,IFMOSTP,MONCOR,NAT1E,
     *                NEEDR,MODRST,NORBPROJ,NUNESP,ISKIPESP,ICURUNT
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
C
      PARAMETER (ZERO=0.0D+00)
C
      DATA RHF,UHF,ROHF,GVB/8HRHF     ,8HUHF     ,8HROHF    ,8HGVB     /
C
C     ----- COMPRESS DENSITY OVER AOS TO DENSITY OVER SHELLS -----
C        PROCESS ONLY SHELLS FROM THE EXTERNAL MONOMER FOR FMO.
C
      IJSH=0
      ALPHA = SCFTYP.EQ.RHF  .OR.  SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
      BETA  =                      SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
      FMOESP=NFG.NE.0.AND.NCURSH.NE.0
      ISH0=1
      JSH0=1
      NEAO=0
      IF(FMOESP) THEN
         ISH0=NCURSH+1
         JSH0=NCURSH+1
         NEAO=L1
      ENDIF
      IJSH=0
      DO 240 ISH=ISH0,NSHELL
         MINI = KLOC(ISH)
         MAXI = MINI + KMAX(ISH) - KMIN(ISH)
         DO 230 JSH=JSH0,ISH
            MINJ = KLOC(JSH)
            MAXJ = MINJ + KMAX(JSH) - KMIN(JSH)
            IJSH = IJSH+1
            DMAX = ZERO
            DO 130 I=MINI,MAXI
               IF(ISH.EQ.JSH) MAXJ=I
               DO 120 J=MINJ,MAXJ
                  IJ = IA(I-NEAO) + J-NEAO
                  IF(ALPHA) THEN
                     IF(ABS(DA(IJ)).GT.DMAX) DMAX = ABS(DA(IJ))
                  END IF
                  IF(BETA) THEN
                     IF(ABS(DB(IJ)).GT.DMAX) DMAX = ABS(DB(IJ))
                  END IF
                  IF(SCFTYP.EQ.GVB) THEN
                     DO 110 IHAM=1,NFO
                        IF(ABS(DGVB(IJ,IHAM)).GT.DMAX)
     *                               DMAX=ABS(DGVB(IJ,IHAM))
  110                CONTINUE
                  END IF
  120          CONTINUE
  130       CONTINUE
            DSH(IJSH) = DMAX
  230    CONTINUE
  240 CONTINUE
      RETURN
      END
C*MODULE SCFLIB  *DECK SOGRAD
      SUBROUTINE SOGRAD(GRAD,F,C,WRK,NPR,NA,L0,L1,ORBGRD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION GRAD(NPR),F(*),C(L1,L0),WRK(L1)
C
      PARAMETER (ZERO=0.0D+00, MXSEQ=150)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,PARR
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
C     -----CALCULATE ORBITAL GRADIENT <OCC|FMO|VIRT> -----
C     WHICH IS C-DAGGER * F * C  TO FORM GRADIENT VECTOR -GRAD-
C          F      - FOCK MATRIX IN AO BASIS (TRIANGULAR)
C          C      - MATRIX OF MO COEFFICIENTS
C          L0     - NUMBER OF LINEARLY INDEPENDENT BASIS FUNCS (MO'S)
C          L1     - NUMBER OF BASIS FUNCTIONS (AO'S)
C          NA     - NUMBER OF OCCUPIED ORBITALS
C          NPR    - NUMBER OF ROTATION PARAMETERS
C          ORBGRD - MAXIMUM ORBITAL GRADIENT COMPONENT
C
      SMALL = 1.0D-03 * CONVHF
      NVIR = L0-NA
      CALL VCLR(WRK,1,L1)
      KG=0
      PARR = GOPARR  .AND.  NA.GT.MXSEQ
C
C          MULTIPLY NA (OCCUPIED) FIRST ROWS OF C-DAGGER BY -F-
C
      DO 300 JO = 1,NA
C
C     ----- GO PARALLEL! -----
C
         IF(PARR) THEN
            IF(MOD(JO,NPROC).NE.ME) GO TO 300
         END IF
C
         IK = 0
         DO 150 I = 1,L1
            IM1 = I-1
            DUM = ZERO
            CDUM = C(I,JO)
            IF (IM1.GT.0) THEN
               DO 100 K = 1,IM1
                  IK = IK+1
                  WRK(K) = WRK(K)+F(IK)*CDUM
                  DUM = DUM+F(IK)*C(K,JO)
  100          CONTINUE
            END IF
            IK = IK+1
            WRK(I) = DUM+F(IK)*CDUM
  150    CONTINUE
C
C           MULTIPLY THESE PRODUCT ROWS BY NVIR COLUMNS IN -C-,
C                   CORRESPONDING ONLY TO VIRTUAL ORBITALS
C
         KG = (JO-1)*NVIR
         DO 200 JV = NA+1,L0
            KG=KG+1
            GRAD(KG) = -DDOT(L1,WRK,1,C(1,JV),1)
            IF (ABS(GRAD(KG)).LT.SMALL) GRAD(KG)=ZERO
  200    CONTINUE
  300 CONTINUE
C
C     ----- GO PARALLEL! -----
C     COLLECT COLUMNS FROM EACH NODE TO MASTER, THEN BROADCAST TO ALL
C
      IF(PARR) THEN
         DO 350 I=1,NA
            IFROM = MOD(I,NPROC)
            IF(IFROM.EQ.MASTER) GO TO 350
            II= (I-1)*NVIR+1
            IF(MASWRK) THEN
               CALL DDI_RECV(GRAD(II),8*NVIR,IFROM)
            ELSE
               IF(IFROM.EQ.ME) CALL DDI_SEND(GRAD(II),8*NVIR,MASTER)
            END IF
  350    CONTINUE
         CALL DDI_BCAST(1011,'F',GRAD,NA*NVIR,MASTER)
      END IF
C
      IMAX   = IDAMAX(NPR,GRAD,1)
      ORBGRD = ABS(GRAD(IMAX))
      RETURN
      END
C*MODULE SCFLIB  *DECK SOGRRO
      SUBROUTINE SOGRRO(GRO,GRA,GRB,NPR,NPRA,NPRB,L1,NA,NB,ORBGRD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION GRO(NPR),GRA(NPRA),GRB(NPRB)
C
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
C     FORM ROHF GRADIENT VECTOR FROM ALPHA AND BETA GRADIENTS
C         L1  = NUMBER OF BASIS FUNCTIONS
C         NA  = NUMBER OF ALPHA MO'S
C         NB  = NUMBER OF BETA MO'S
C         NPR = NUMBER OF ROHF PAIR ROTATIONS
C
      SMALL = 1.0D-03 * CONVHF
      K=0
      KA=0
      KB=0
      DO 110 J=1,NA
         IF (J.LE.NB) THEN
            N=NB+1
         ELSE
            N=NA+1
         END IF
         DO 100 I=N,L1
            K=K+1
            IF (J.LE.NB) THEN
               KB=KB+1
               IF (I.LE.NA) THEN
                  GRO(K)=GRB(KB)
               ELSE
                  KA=KA+1
                  GRO(K) = (GRB(KB)+GRA(KA))*0.5D+00
                  IF (ABS(GRO(K)).LT.SMALL) GRO(K)=0.0D+00
               END IF
            ELSE
               KA=KA+1
               GRO(K) = GRA(KA)
            END IF
  100    CONTINUE
  110 CONTINUE
C
      IMAX   = IDAMAX(NPR,GRO,1)
      ORBGRD = ABS(GRO(IMAX))
C
      RETURN
      END
C*MODULE SCFLIB  *DECK SOHESS
      SUBROUTINE SOHESS(HSTART,EIG,NPR,L1,NA,NB)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION HSTART(NPR),EIG(L1)
C
      PARAMETER (ONE=1.0D+00, PT05=0.5D-02, TWENTY=ONE/PT05)
C
C     ----- SET APPROXIMATE INVERSE HESSIAN -----
C     DIAGONAL ELEMENTS OF INITIAL ORBITAL HESSIAN -HSTART-
C     ARE SET USING ONLY THE ENERGY EIGENVALUES -EIG-
C
C         L1  = NUMBER OF BASIS FUNCTIONS
C         NA  = NUMBER OF ALPHA MO'S
C         NB  = NUMBER OF BETA MO'S
C         NPR = NUMBER OF PAIR ROTATIONS
C
      K=0
      DO 110 J=1,NA
         IF (J.LE.NB) THEN
            N=NB+1
         ELSE
            N=NA+1
         END IF
         DO 100 I=N,L1
            K=K+1
            DE = EIG(I)-EIG(J)
            IF(DE.GT.PT05) THEN
               HSTART(K)=ONE/DE
            ELSE
               HSTART(K)=TWENTY
            END IF
  100    CONTINUE
  110 CONTINUE
      RETURN
      END
C*MODULE SCFLIB  *DECK SONEWT
      SUBROUTINE SONEWT(HSTART,GRAD,PGRAD,DISPLI,DGRAD,DISPL,UPDT,
     *                  DISPLN,DGRADI,UPDTI,ORBGRD,NPR,ITSO,NFT15)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION HSTART(NPR),GRAD(NPR),PGRAD(NPR),DISPLI(NPR),
     *          DGRAD(NPR),DISPL(NPR),UPDT(NPR),DISPLN(NPR),
     *          DGRADI(NPR),UPDTI(NPR)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SVDSKW,GPSAVE
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
      PARAMETER (ONE=1.0D+00)
C
C     ----- THIS ROUTINE PERFORMS INVERSE HESSIAN MATRIX UPDATES -----
C                USING THE DIRECT UPDATE PROCEDURE GIVEN IN
C           T.H.FISCHER, J.ALMLOF.  J.PHYS.CHEM.  1992,96,9768-74
C
C        HSTART - APPROXIMATE INITIAL DIAGONAL INVERSE HESSIAN
C        GRAD   - GRADIENT VECTOR AT CURRENT POINT
C        PGRAD  - GRADIENT FROM PREVIOUS POINT
C        DISPLI - ON EXIT, THE DISPLACEMENT VECTOR
C        ALL OTHER ARGUMENTS ARE SCRATCH STORAGE.
C
      SMALL = 1.0D-03 * CONVHF
      TOOBIG = 1.0D+00
      BIGROT = 0.1D+00
C
C     THIS ROUTINE CONTAINS WORK PROPORTIONAL TO N**2, AND SOME I/O,
C     AND THUS IS MOST EFFICIENTLY DONE BY THE MASTER NODE ONLY.
C
      SVDSKW = DSKWRK
      GPSAVE = GOPARR
      DSKWRK = .FALSE.
      GOPARR = .FALSE.
      IF(.NOT.MASWRK) GO TO 500
C
C     INITIALIZE DISPLACEMENT VECTOR
C
      DO 100 I=1,NPR
         DISPLN(I)=HSTART(I)*GRAD(I)
  100 CONTINUE
C
      IF (ITSO.EQ.1) GO TO 400
C
      IF (ORBGRD.LT.SMALL) THEN
         ITSO = ITSO - 1
         GO TO 400
      END IF
C
      DO 120 I=1,NPR
         DGRADI(I)= GRAD(I)-PGRAD(I)
         UPDTI(I) = HSTART(I)*DGRADI(I)
  120 CONTINUE
C
      IF (ITSO.EQ.2) GO TO 300
C
C     CALCULATE VECTOR RECURSIVELY
C
      CALL SEQREW(NFT15)
C
      DO 200 K=1,ITSO-2
C
         CALL SQREAD(NFT15,DISPL,NPR)
         CALL SQREAD(NFT15,DGRAD,NPR)
         CALL SQREAD(NFT15, UPDT,NPR)
C
         S1 = DDOT(NPR, DISPL,1, DGRAD,1)
         S2 = DDOT(NPR, DGRAD,1,  UPDT,1)
         S3 = DDOT(NPR, DISPL,1,  GRAD,1)
         S4 = DDOT(NPR,  UPDT,1,  GRAD,1)
         S5 = DDOT(NPR, DISPL,1,DGRADI,1)
         S6 = DDOT(NPR,  UPDT,1,DGRADI,1)
C
         S1=ONE/S1
         S2=ONE/S2
         T=ONE+S1/S2
         T2=S1*S3
         T4=S1*S5
         T1=T*T2-S1*S4
         T3=T*T4-S1*S6
C
         CALL DAXPY(NPR, T1,DISPL,1,DISPLN,1)
         CALL DAXPY(NPR,-T2, UPDT,1,DISPLN,1)
         CALL DAXPY(NPR, T3,DISPL,1, UPDTI,1)
         CALL DAXPY(NPR,-T4, UPDT,1, UPDTI,1)
C
  200 CONTINUE
C
C     UPDATE VECTOR
C
  300 CONTINUE
      S1 = DDOT(NPR,DISPLI,1,DGRADI,1)
      S2 = DDOT(NPR,DGRADI,1, UPDTI,1)
      S3 = DDOT(NPR,DISPLI,1,  GRAD,1)
      S4 = DDOT(NPR, UPDTI,1,  GRAD,1)
C
      S1=ONE/S1
      S2=ONE/S2
      T=ONE+S1/S2
      T2=S1*S3
      T1=T*T2-S1*S4
C
      CALL DAXPY(NPR, T1,DISPLI,1,DISPLN,1)
      CALL DAXPY(NPR,-T2, UPDTI,1,DISPLN,1)
C
      CALL SQWRIT(NFT15,DISPLI,NPR)
      CALL SQWRIT(NFT15,DGRADI,NPR)
      CALL SQWRIT(NFT15, UPDTI,NPR)
C
  400 CONTINUE
      DO 410 I=1,NPR
         DISPLI(I) = -DISPLN(I)
  410 CONTINUE
C
C        SCALE DISPLACEMENT SO THAT SQCDF DOESN'T EXCEED 0.1
C
      SQCDF = SQRT(DDOT(NPR,DISPLI,1,DISPLI,1)/NPR)
      IF(SQCDF.GT.TOOBIG  .AND.  ITSO.GT.5) THEN
         IF(MASWRK) WRITE(IW,9010) SQCDF
         ISTAT=1
         CALL DDI_BCAST(1016,'I',ISTAT,1,MASTER)
         IF(NFG.EQ.0) THEN
           CALL ABRT
         ELSE
C          ABORTING HERE RESULTS IN DEADLOCKS FOR GDDI
           ORBGRD=-123
           DSKWRK = SVDSKW
           GOPARR = GPSAVE
           RETURN
         ENDIF
      END IF
      IF (SQCDF.GT.BIGROT) THEN
         IF(MASWRK) WRITE(IW,9020) SQCDF
         SCAL=SQRT(BIGROT/SQCDF)
         CALL DSCAL(NPR,SCAL,DISPLI,1)
      END IF
      ISTAT=0
      CALL DDI_BCAST(1016,'I',ISTAT,1,MASTER)
C
C        GIVE DISPLACEMENT TO ANY OTHER NODES
C
  500 CONTINUE
      DSKWRK = SVDSKW
      GOPARR = GPSAVE
C
C       OTHER NODES MUST CHECK TO SEE IF THE SOSCF RAN OK ON MASTER
C
      IF(.NOT.MASWRK) CALL DDI_BCAST(1016,'I',ISTAT,1,MASTER)
      IF(ISTAT.EQ.1) THEN
         IF(NFG.EQ.0) THEN
           CALL ABRT
         ELSE
           ORBGRD=-123
           RETURN
         ENDIF
      END IF
C
      IF(GOPARR) THEN
         CALL DDI_BCAST(1015,'F',DISPLI,NPR,MASTER)
      END IF
      RETURN
C
 9010 FORMAT(1X,'SOSCF ENCOUNTERS A SERIOUS PROBLEM IN -SONEWT-'/
     *       1X,'THE ROTATION ANGLE VECTOR HAS A HUGE NORM, SQCDF=',
     *          1P,E12.3/1X,'REEXAMINE STARTING VECTORS, ',
     *             'APPROPRIATENESS OF YOUR SCFTYP, ETC...')
 9020 FORMAT(1X,'SOSCF IS SCALING ROTATION ANGLE MATRIX, SQCDF=',F12.6)
      END
C*MODULE SCFLIB  *DECK SOTRAN
      SUBROUTINE SOTRAN(X,C,G,WRK,NPR,L0,L1,NA,NB,ORBGRD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION X(NPR),C(L1,L0),G(L0,L0),WRK(L0)
C
      PARAMETER (ZERO=0.0D+00, HALF=0.5D+00, ONE=1.0D+00, MXSEQ=150)
      LOGICAL GOPARR,DSKWRK,MASWRK,PARR
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- APPROXIMATE SECOND ORDER ORBITAL TRANSFORMATION -----
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          NA     = NUMBER OF ALPHA MO'S
C          NB     = NUMBER OF BETA MO'S
C          NPR    = NUMBER OF ROTATION PARAMETERS ( NPR=NOCC*NVIR)
C
C     ORIGINAL ROUTINE WRITTEN BY GALINA CHABAN AT ISU, 1995
C     PARALLELIZATION IS BY KAZUYA ISHIMURA AT IMS, 2004
C
      CALL VCLR(G,1,L0*L0)
      PARR = GOPARR  .AND.  L0.GT.MXSEQ
C
C     ----- FORM MATRIX G = I + X + ( X*X/2 ) -----
C
      K=0
      DO 110 J=1,NA
         IF (J.LE.NB) THEN
            N=NB+1
         ELSE
            N=NA+1
         END IF
         DO 100 I=N,L0
            K=K+1
            G(J,I)= X(K)
            G(I,J)=-X(K)
  100    CONTINUE
  110 CONTINUE
C
C        TESTS SHOWED THE FULL SECOND ORDER EXPANSION IS OF LITTLE USE
C        SO WE ALWAYS SKIP IT, BUT SAVE THE CODE FOR FUTURE REFERENCE.
C        (SECOND ORDER EXPANSION CODE WORKS ONLY FOR NA=NB)
C
      IF(NA.NE.NB) GO TO 300
      IF(ORBGRD.GT.1.0D-20) GO TO 300
C
      NOCC=NA
      NVIR=L0-NOCC
C
      DO 230 I=1,NOCC
         DO 220 J=I,NOCC
            DUM=ZERO
            DO 210 K=1,NVIR
               LI=NVIR*(I-1)+K
               LJ=NVIR*(J-1)+K
               DUM=DUM+X(LI)*X(LJ)
  210       CONTINUE
            DUM = -DUM*HALF
            G(I,J)= DUM
            G(J,I)= DUM
  220    CONTINUE
  230 CONTINUE
C
      DO 260 I=1,NVIR
         DO 250 J=I,NVIR
            DUM=ZERO
            DO 240 K=1,NOCC
               LI=NVIR*(K-1)+I
               LJ=LI-I+J
               DUM=DUM+X(LI)*X(LJ)
  240       CONTINUE
            DUM = -DUM*HALF
            I1=NOCC+I
            J1=NOCC+J
            G(I1,J1)= DUM
            G(J1,I1)= DUM
  250    CONTINUE
  260 CONTINUE
C
  300 CONTINUE
      DO 310 I=1,L0
         G(I,I) = G(I,I) + ONE
  310 CONTINUE
C
C     ----- ORTHONORMALIZE THE TRANSFORMATION MATRIX -----
C
      IF(PARR) THEN
         DO 520 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 520
            IP1=I+1
            DO 510 J=IP1,L0
               IF(MOD(J,NPROC).NE.ME) GO TO 510
               DUM = DDOT(L0,G(1,I),1,G(1,J),1)
               CALL DAXPY(L0,-DUM,G(1,I),1,G(1,J),1)
  510       CONTINUE
            CALL DDI_BCAST(1011,'F',G(1,I+1),L0,MOD(I+1,NPROC))
  520    CONTINUE
      ELSE
         DO 540 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 540
            IP1=I+1
            DO 530 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)
  530       CONTINUE
  540    CONTINUE
      END IF
C
C     ----- ROTATE THE ORBITALS -C- BY TRANSFORMATION -G- -----
C
      IF(PARR) THEN
         DO 650 I=1,L1
            DO 610 K=1,L0
               WRK(K)=C(I,K)
  610       CONTINUE
            DO 630 J=1,L0
               IF(MOD(J,NPROC).NE.ME) GO TO 630
               DUM=ZERO
               DO 620 K=1,L0
                  DUM=DUM+WRK(K)*G(K,J)
  620          CONTINUE
               C(I,J)=DUM
  630       CONTINUE
  650    CONTINUE
C
         DO 660 J=1,L0
            JFROM = MOD(J,NPROC)
            IF(JFROM.EQ.MASTER) GO TO 660
            IF(MASWRK) THEN
               CALL DDI_RECV(C(1,J),8*L1,JFROM)
            ELSE
               IF(JFROM.EQ.ME) CALL DDI_SEND(C(1,J),8*L1,MASTER)
            END IF
  660    CONTINUE
         CALL DDI_BCAST(1012,'F',C,L1*L0,MASTER)
      ELSE
         DO 750 I=1,L1
            DO 710 K=1,L0
               WRK(K)=C(I,K)
  710       CONTINUE
            DO 730 J=1,L0
               DUM=ZERO
               DO 720 K=1,L0
                  DUM=DUM+WRK(K)*G(K,J)
  720          CONTINUE
               C(I,J)=DUM
  730       CONTINUE
  750    CONTINUE
      END IF
      RETURN
      END
C*MODULE SCFLIB  *DECK SUBMT
      SUBROUTINE SUBMT(A,LDA)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(LDA,LDA)
      PARAMETER (ZERO=0.0D+00)
C
C        SET MATRIX A = A - A-TRANSPOSE
C
      DO 120 J=1,LDA
         JM1=J-1
         DO 110 I=1,JM1
            DIFF   = A(I,J)-A(J,I)
            A(I,J) = DIFF
            A(J,I) =-DIFF
  110    CONTINUE
         A(J,J)=ZERO
  120 CONTINUE
      RETURN
      END
C*MODULE SCFLIB  *DECK SYMH
      SUBROUTINE SYMH(F,H,IA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL IANDJ
C
      DIMENSION F(*),H(*),IA(*)
      DIMENSION MI(48)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
      COMMON /HSYM  / T(35,35),MINI,MAXI,LIT,MINJ,MAXJ,LJT,NTR
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TT(432),INVT(48),NT
C
      PARAMETER (ONE=1.0D+00)
C
C     ----- SYMMETRIZE THE SKELETON FOCK MATRIX
C
      IF (NT .EQ. 1) RETURN
C
      L2 = (NUM*NUM+NUM)/2
      CALL VCLR(H,1,L2)
C
C     ----- FIND A BLOCK (I,J)
C
      DO 520 II = 1,NSHELL
      DO 140 ITR = 1,NT
         ISH = MAPSHL(II,ITR)
         IF (ISH .GT. II) GO TO 520
         MI(ITR) = ISH
  140 CONTINUE
      LIT = KTYPE(II)
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II)-MINI
      DO 500 JJ = 1,II
      DO 200 ITR = 1,NT
      JSH = MAPSHL(JJ,ITR)
      IF (JSH .GT. II) GO TO 500
      ISH = MI(ITR)
      IF (ISH .GE. JSH) GO TO 180
      N = ISH
      ISH = JSH
      JSH = N
  180 IF (ISH .EQ. II .AND. JSH .GT. JJ) GO TO 500
  200 CONTINUE
      LJT = KTYPE(JJ)
      MINJ = KMIN(JJ)
      MAXJ = KMAX(JJ)
      LOCJ = KLOC(JJ)-MINJ
      IANDJ = II .EQ. JJ
      JMAX = MAXJ
C
C     ----- FIND THE EQUIVALENT BLOCKS -----
C     ----- TRANSFER EQUIVALENT BLOCK INTO T-MATRIX
C     ----- COMPUTE (R) T (R)
C     ----- PUT THE RESULT BACK INTO THE (I,J) BLOCK OF THE H-MATRIX
C
      DO 300 ITR = 1,NT
      NTR = ITR
      KK = MI(ITR)
      LL = MAPSHL(JJ,ITR)
      LOCK = KLOC(KK)-KMIN(KK)
      LOCL = KLOC(LL)-KMIN(LL)
      DO 260 K = MINI,MAXI
      LCK = LOCK+K
      IF (IANDJ) JMAX = K
      DO 260 L = MINJ,JMAX
      IF (LL .GT. KK) GO TO 220
      KL = IA(LCK)+LOCL+L
      GO TO 240
C
  220 KL = IA(LOCL+L)+LCK
  240 T(K,L) = F(KL)
      IF (IANDJ) T(L,K) = F(KL)
  260 CONTINUE
      CALL RHR
      DO 280 I = MINI,MAXI
      LCI = IA(LOCI+I)+LOCJ
      IF (IANDJ) JMAX = I
      DO 280 J = MINJ,JMAX
      IJ = LCI+J
  280 H(IJ) = H(IJ)+T(I,J)
  300 CONTINUE
C
C     ----- FOR EACH BLOCK (K,L) EQUIVALENT TO (I,J)
C     ----- FIND THE TRANSFORMATION THAT MAPS (K,L) INTO (I,J)
C     ----- COMPUTE (R) T (R)
C     ----- PUT THE RESULT BACK INTO THE (K,L) BLOCK OF THE H-MATRIX
C
      DO 480 ITR = 1,NT
      IF (ITR .EQ. 1) GO TO 480
      KK = MI(ITR)
      LL = MAPSHL(JJ,ITR)
      IF (KK .GE. LL) GO TO 320
      K = LL
      L = KK
      GO TO 340
C
  320 K = KK
      L = LL
  340 IF (K .EQ. II .AND. L .EQ. JJ) GO TO 480
      NTR = ITR+1
      IF (NTR .GT. NT) GO TO 400
      DO 380 IT = NTR,NT
      I = MI(IT)
      J = MAPSHL(JJ,IT)
      IF (I .GE. J) GO TO 360
      IJ = I
      I = J
      J = IJ
  360 IF (I .EQ. K .AND. J .EQ. L) GO TO 480
  380 CONTINUE
  400 CONTINUE
      NTR = INVT(ITR)
      DO 420 I = MINI,MAXI
      LCI = IA(LOCI+I)+LOCJ
      IF (IANDJ) JMAX = I
      DO 420 J = MINJ,JMAX
      T(I,J) = H(LCI+J)
      IF (IANDJ) T(J,I) = H(LCI+J)
  420 CONTINUE
      CALL RHR
      LOCK = KLOC(KK)-KMIN(KK)
      LOCL = KLOC(LL)-KMIN(LL)
      DO 460 K = MINI,MAXI
      LCK = LOCK+K
      IF (IANDJ) JMAX = K
      DO 460 L = MINJ,JMAX
      IF (LL .GT. KK) GO TO 440
      KL = IA(LCK)+LOCL+L
      GO TO 460
C
  440 KL = IA(LOCL+L)+LCK
  460 H(KL) = T(K,L)
  480 CONTINUE
  500 CONTINUE
  520 CONTINUE
      DUM = ONE/NT
      DO 540 I = 1,L2
  540 F(I) = H(I)*DUM
      RETURN
      END
C*MODULE SCFLIB  *DECK TRMAT
      SUBROUTINE TRMAT
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION PNAME(3),DNAME(6),FNAME(10),GNAME(15)
C
      PARAMETER (MXSH=5000, MXATM=2000)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TT(432),INVT(48),NT
      COMMON /TRANSF/ XSMAL,YSMAL,ZSMAL,XNEW,YNEW,ZNEW,XP,YP,ZP
      COMMON /FRAME / U1,U2,U3,V1,V2,V3,W1,W2,W3,X0,Y0,Z0
C
      PARAMETER (ONE=1.0D+00, SQRT3=1.73205080756888D+00,
     *           SQRT5=2.23606797749979D+00,
     *           SQRT7=2.64575131106459D+00)
C
      DATA PNAME /4HX   ,4HY   ,4HZ   /
      DATA DNAME /4HXX  ,4HYY  ,4HZZ  ,4HXY  ,4HXZ  ,4HYZ  /
      DATA FNAME /4HXXX ,4HYYY ,4HZZZ ,4HXXY ,4HXXZ ,4HYYX ,
     *            4HYYZ ,4HZZX ,4HZZY ,4HXYZ /
      DATA GNAME /4HXXXX,4HYYYY,4HZZZZ,4HXXXY,4HXXXZ,4HYYYX,
     *            4HYYYZ,4HZZZX,4HZZZY,4HXXYY,4HXXZZ,4HYYZZ,
     *            4HXXYZ,4HYYXZ,4HZZXY/
C
C     ----- CALCULATE TRANSFORMS OF P, D, F, AND G FUNCTIONS
C           FOR ALL SYMETRY OPERATIONS.
C
      XS = ONE
      YS = ONE
      ZS = ONE
      X = X0+ONE
      Y = Y0
      Z = Z0
      CALL LOCAL(X,Y,Z,XS,YS,ZS)
      XSMAL = XS
      YSMAL = YS
      ZSMAL = ZS
      DO 100 IT = 1,NT
      NN = 9*(IT-1)
      CALL TRANS(NN)
      CALL ROT
      N = 3*(IT-1)
      PTR(1,N+1) = XP-X0
      PTR(2,N+1) = YP-Y0
      PTR(3,N+1) = ZP-Z0
  100 CONTINUE
      X = X0
      Y = Y0+ONE
      Z = Z0
      CALL LOCAL(X,Y,Z,XS,YS,ZS)
      XSMAL = XS
      YSMAL = YS
      ZSMAL = ZS
      DO 120 IT = 1,NT
      NN = 9*(IT-1)
      CALL TRANS(NN)
      CALL ROT
      N = 3*(IT-1)
      PTR(1,N+2) = XP-X0
      PTR(2,N+2) = YP-Y0
      PTR(3,N+2) = ZP-Z0
  120 CONTINUE
      X = X0
      Y = Y0
      Z = Z0+ONE
      CALL LOCAL(X,Y,Z,XS,YS,ZS)
      XSMAL = XS
      YSMAL = YS
      ZSMAL = ZS
      DO 140 IT = 1,NT
      NN = 9*(IT-1)
      CALL TRANS(NN)
      CALL ROT
      N = 3*(IT-1)
      PTR(1,N+3) = XP-X0
      PTR(2,N+3) = YP-Y0
      PTR(3,N+3) = ZP-Z0
  140 CONTINUE
      DO 640 IT = 1,NT
      NP = 3*(IT-1)
      ND = 6*(IT-1)
      NF = 10*(IT-1)
      NG = 15*(IT-1)
      DO 300 I = 1,6
      GO TO (160,180,200,220,240,260) ,I
  160 J = 1
      K = 1
      GO TO 280
C
  180 J = 2
      K = 2
      GO TO 280
C
  200 J = 3
      K = 3
      GO TO 280
C
  220 J = 1
      K = 2
      GO TO 280
C
  240 J = 1
      K = 3
      GO TO 280
C
  260 J = 2
      K = 3
  280 DTR(1,ND+I) = PTR(1,NP+J)*PTR(1,NP+K)
      DTR(2,ND+I) = PTR(2,NP+J)*PTR(2,NP+K)
      DTR(3,ND+I) = PTR(3,NP+J)*PTR(3,NP+K)
      DTR(4,ND+I) = PTR(1,NP+J)*PTR(2,NP+K) +PTR(2,NP+J)*PTR(1,NP+K)
      DTR(5,ND+I) = PTR(1,NP+J)*PTR(3,NP+K) +PTR(3,NP+J)*PTR(1,NP+K)
      DTR(6,ND+I) = PTR(2,NP+J)*PTR(3,NP+K) +PTR(3,NP+J)*PTR(2,NP+K)
  300 CONTINUE
      DO 440 I=1,10
      GO TO (320,330,340,350,360,370,380,390,400,410),I
  320 J=1
      K=1
      GO TO 420
  330 J=2
      K=2
      GO TO 420
  340 J=3
      K=3
      GO TO 420
  350 J=1
      K=2
      GO TO 420
  360 J=1
      K=3
      GO TO 420
  370 J=2
      K=1
      GO TO 420
  380 J=2
      K=3
      GO TO 420
  390 J=3
      K=1
      GO TO 420
  400 J=3
      K=2
      GO TO 420
  410 J=4
      K=3
  420 FTR(1,NF+I)=DTR(1,ND+J)*PTR(1,NP+K)
      FTR(2,NF+I)=DTR(2,ND+J)*PTR(2,NP+K)
      FTR(3,NF+I)=DTR(3,ND+J)*PTR(3,NP+K)
      FTR(4,NF+I)=DTR(1,ND+J)*PTR(2,NP+K)
     1           +DTR(4,ND+J)*PTR(1,NP+K)
      FTR(5,NF+I)=DTR(1,ND+J)*PTR(3,NP+K)
     1           +DTR(5,ND+J)*PTR(1,NP+K)
      FTR(6,NF+I)=DTR(2,ND+J)*PTR(1,NP+K)
     1           +DTR(4,ND+J)*PTR(2,NP+K)
      FTR(7,NF+I)=DTR(2,ND+J)*PTR(3,NP+K)
     1           +DTR(6,ND+J)*PTR(2,NP+K)
      FTR(8,NF+I)=DTR(3,ND+J)*PTR(1,NP+K)
     1           +DTR(5,ND+J)*PTR(3,NP+K)
      FTR(9,NF+I)=DTR(3,ND+J)*PTR(2,NP+K)
     1           +DTR(6,ND+J)*PTR(3,NP+K)
      FTR(10,NF+I)=DTR(4,ND+J)*PTR(3,NP+K)
     1           + DTR(5,ND+J)*PTR(2,NP+K)
     2           + DTR(6,ND+J)*PTR(1,NP+K)
  440 CONTINUE
      DO 620 I=1,15
      GO TO (460,470,480,490,500,510,520,530,
     1       540,550,560,570,580,590,600),I
  460 J=1
      K=1
      GO TO 610
  470 J=2
      K=2
      GO TO 610
  480 J=3
      K=3
      GO TO 610
  490 J=1
      K=2
      GO TO 610
  500 J=1
      K=3
      GO TO 610
  510 J=2
      K=1
      GO TO 610
  520 J=2
      K=3
      GO TO 610
  530 J=3
      K=1
      GO TO 610
  540 J=3
      K=2
      GO TO 610
  550 J=4
      K=2
      GO TO 610
  560 J=5
      K=3
      GO TO 610
  570 J=7
      K=3
      GO TO 610
  580 J=4
      K=3
      GO TO 610
  590 J=6
      K=3
      GO TO 610
  600 J=8
      K=2
  610 GTR(1,NG+I)=FTR(1,NF+J)*PTR(1,NP+K)
      GTR(2,NG+I)=FTR(2,NF+J)*PTR(2,NP+K)
      GTR(3,NG+I)=FTR(3,NF+J)*PTR(3,NP+K)
      GTR(4,NG+I)=FTR(1,NF+J)*PTR(2,NP+K)
     1           +FTR(4,NF+J)*PTR(1,NP+K)
      GTR(5,NG+I)=FTR(1,NF+J)*PTR(3,NP+K)
     1           +FTR(5,NF+J)*PTR(1,NP+K)
      GTR(6,NG+I)=FTR(2,NF+J)*PTR(1,NP+K)
     1           +FTR(6,NF+J)*PTR(2,NP+K)
      GTR(7,NG+I)=FTR(2,NF+J)*PTR(3,NP+K)
     1           +FTR(7,NF+J)*PTR(2,NP+K)
      GTR(8,NG+I)=FTR(3,NF+J)*PTR(1,NP+K)
     1           +FTR(8,NF+J)*PTR(3,NP+K)
      GTR(9,NG+I)=FTR(3,NF+J)*PTR(2,NP+K)
     1           +FTR(9,NF+J)*PTR(3,NP+K)
      GTR(10,NG+I)=FTR(4,NF+J)*PTR(2,NP+K)
     1           + FTR(6,NF+J)*PTR(1,NP+K)
      GTR(11,NG+I)=FTR(5,NF+J)*PTR(3,NP+K)
     1           + FTR(8,NF+J)*PTR(1,NP+K)
      GTR(12,NG+I)=FTR(7,NF+J)*PTR(3,NP+K)
     1           + FTR(9,NF+J)*PTR(2,NP+K)
      GTR(13,NG+I)=FTR(4,NF+J)*PTR(3,NP+K)
     1           + FTR(5,NF+J)*PTR(2,NP+K)
     2           +FTR(10,NF+J)*PTR(1,NP+K)
      GTR(14,NG+I)=FTR(6,NF+J)*PTR(3,NP+K)
     1           + FTR(7,NF+J)*PTR(1,NP+K)
     2           +FTR(10,NF+J)*PTR(2,NP+K)
      GTR(15,NG+I)=FTR(8,NF+J)*PTR(2,NP+K)
     1           + FTR(9,NF+J)*PTR(1,NP+K)
     2           +FTR(10,NF+J)*PTR(3,NP+K)
  620 CONTINUE
  640 CONTINUE
      IF (NORMF .EQ. 1 .AND. NORMP .EQ. 1) GO TO 750
      DO 740 IT = 1,NT
      ND = 6*(IT-1)
      NF = 10*(IT-1)
      NG = 15*(IT-1)
      DO 660 I = 1,6
      IF (I .GT. 3) GO TO 650
      DTR(4,ND+I) = DTR(4,ND+I)/SQRT3
      DTR(5,ND+I) = DTR(5,ND+I)/SQRT3
      DTR(6,ND+I) = DTR(6,ND+I)/SQRT3
      GO TO 660
C
  650 DTR(1,ND+I) = DTR(1,ND+I)*SQRT3
      DTR(2,ND+I) = DTR(2,ND+I)*SQRT3
      DTR(3,ND+I) = DTR(3,ND+I)*SQRT3
  660 CONTINUE
      DO 690 I=1,10
      IF(I.GT.3) GO TO 670
      FTR(4,NF+I)=FTR(4,NF+I)/SQRT5
      FTR(5,NF+I)=FTR(5,NF+I)/SQRT5
      FTR(6,NF+I)=FTR(6,NF+I)/SQRT5
      FTR(7,NF+I)=FTR(7,NF+I)/SQRT5
      FTR(8,NF+I)=FTR(8,NF+I)/SQRT5
      FTR(9,NF+I)=FTR(9,NF+I)/SQRT5
      FTR(10,NF+I)=FTR(10,NF+I)/(SQRT5*SQRT3)
      GO TO 690
  670 IF(I.GT.9) GO TO 680
      FTR(1,NF+I)=FTR(1,NF+I)*SQRT5
      FTR(2,NF+I)=FTR(2,NF+I)*SQRT5
      FTR(3,NF+I)=FTR(3,NF+I)*SQRT5
      FTR(10,NF+I)=FTR(10,NF+I)/SQRT3
      GO TO 690
  680 FTR(1,NF+I)=FTR(1,NF+I)*SQRT5*SQRT3
      FTR(2,NF+I)=FTR(2,NF+I)*SQRT5*SQRT3
      FTR(3,NF+I)=FTR(3,NF+I)*SQRT5*SQRT3
      FTR(4,NF+I)=FTR(4,NF+I)*SQRT3
      FTR(5,NF+I)=FTR(5,NF+I)*SQRT3
      FTR(6,NF+I)=FTR(6,NF+I)*SQRT3
      FTR(7,NF+I)=FTR(7,NF+I)*SQRT3
      FTR(8,NF+I)=FTR(8,NF+I)*SQRT3
      FTR(9,NF+I)=FTR(9,NF+I)*SQRT3
  690 CONTINUE
      DO 730 I=1,15
      IF(I.GT.3) GO TO 700
      GTR(4,NG+I)=GTR(4,NG+I)/SQRT7
      GTR(5,NG+I)=GTR(5,NG+I)/SQRT7
      GTR(6,NG+I)=GTR(6,NG+I)/SQRT7
      GTR(7,NG+I)=GTR(7,NG+I)/SQRT7
      GTR(8,NG+I)=GTR(8,NG+I)/SQRT7
      GTR(9,NG+I)=GTR(9,NG+I)/SQRT7
      GTR(10,NG+I)=GTR(10,NG+I)*SQRT3/(SQRT5*SQRT7)
      GTR(11,NG+I)=GTR(11,NG+I)*SQRT3/(SQRT5*SQRT7)
      GTR(12,NG+I)=GTR(12,NG+I)*SQRT3/(SQRT5*SQRT7)
      GTR(13,NG+I)=GTR(13,NG+I)/(SQRT5*SQRT7)
      GTR(14,NG+I)=GTR(14,NG+I)/(SQRT5*SQRT7)
      GTR(15,NG+I)=GTR(15,NG+I)/(SQRT5*SQRT7)
      GO TO 730
  700 IF(I.GT.9) GO TO 710
      GTR(1,NG+I)=GTR(1,NG+I)*SQRT7
      GTR(2,NG+I)=GTR(2,NG+I)*SQRT7
      GTR(3,NG+I)=GTR(3,NG+I)*SQRT7
      GTR(10,NG+I)=GTR(10,NG+I)*SQRT3/SQRT5
      GTR(11,NG+I)=GTR(11,NG+I)*SQRT3/SQRT5
      GTR(12,NG+I)=GTR(12,NG+I)*SQRT3/SQRT5
      GTR(13,NG+I)=GTR(13,NG+I)/SQRT5
      GTR(14,NG+I)=GTR(14,NG+I)/SQRT5
      GTR(15,NG+I)=GTR(15,NG+I)/SQRT5
      GO TO 730
  710 IF(I.GT.12) GO TO 720
      GTR(1,NG+I)=GTR(1,NG+I)*SQRT7*SQRT5/SQRT3
      GTR(2,NG+I)=GTR(2,NG+I)*SQRT7*SQRT5/SQRT3
      GTR(3,NG+I)=GTR(3,NG+I)*SQRT7*SQRT5/SQRT3
      GTR(4,NG+I)=GTR(4,NG+I)*SQRT5/SQRT3
      GTR(5,NG+I)=GTR(5,NG+I)*SQRT5/SQRT3
      GTR(6,NG+I)=GTR(6,NG+I)*SQRT5/SQRT3
      GTR(7,NG+I)=GTR(7,NG+I)*SQRT5/SQRT3
      GTR(8,NG+I)=GTR(8,NG+I)*SQRT5/SQRT3
      GTR(9,NG+I)=GTR(9,NG+I)*SQRT5/SQRT3
      GTR(13,NG+I)=GTR(13,NG+I)/SQRT3
      GTR(14,NG+I)=GTR(14,NG+I)/SQRT3
      GTR(15,NG+I)=GTR(15,NG+I)/SQRT3
      GO TO 730
  720 GTR(1,NG+I)=GTR(1,NG+I)*SQRT7*SQRT5
      GTR(2,NG+I)=GTR(2,NG+I)*SQRT7*SQRT5
      GTR(3,NG+I)=GTR(3,NG+I)*SQRT7*SQRT5
      GTR(4,NG+I)=GTR(4,NG+I)*SQRT5
      GTR(5,NG+I)=GTR(5,NG+I)*SQRT5
      GTR(6,NG+I)=GTR(6,NG+I)*SQRT5
      GTR(7,NG+I)=GTR(7,NG+I)*SQRT5
      GTR(8,NG+I)=GTR(8,NG+I)*SQRT5
      GTR(9,NG+I)=GTR(9,NG+I)*SQRT5
      GTR(10,NG+I)=GTR(10,NG+I)*SQRT3
      GTR(11,NG+I)=GTR(11,NG+I)*SQRT3
      GTR(12,NG+I)=GTR(12,NG+I)*SQRT3
  730 CONTINUE
  740 CONTINUE
  750 CONTINUE
C
C     ----- PRINT MATRICES IF NPRINT.EQ.5 -----
C
      IF (NPRINT .NE. 5) GO TO 860
      IF (MASWRK) THEN
      WRITE (IW,9028)
      DO 840 IT = 1,NT
      WRITE (IW,9088)
      WRITE (IW,9108) IT
      NP = 3*(IT-1)
      WRITE (IW,9048) (PNAME(J),J = 1,3)
      WRITE (IW,9128)
      DO 760 I = 1,3
  760 WRITE (IW,9068) PNAME(I),(PTR(I,NP+J),J = 1,3)
      WRITE (IW,9008)
      ND = 6*(IT-1)
      WRITE (IW,9048) (DNAME(J),J = 1,6)
      WRITE (IW,9128)
      DO 780 I = 1,6
  780 WRITE (IW,9068) DNAME(I),(DTR(I,ND+J),J = 1,6)
      WRITE (IW,9008)
      NF = 10*(IT-1)
      WRITE (IW,9048) (FNAME(J),J = 1,10)
      WRITE (IW,9128)
      DO 800 I = 1,10
  800 WRITE (IW,9068) FNAME(I),(FTR(I,ND+J),J = 1,10)
      WRITE (IW,9008)
      NG = 15*(IT-1)
      WRITE (IW,9048) (GNAME(J),J = 1,15)
      WRITE (IW,9128)
      DO 820 I = 1,15
  820 WRITE (IW,9068) GNAME(I),(GTR(I,ND+J),J = 1,15)
      WRITE (IW,9008)
  840 CONTINUE
      END IF
  860 CONTINUE
      RETURN
C
 9008 FORMAT(//)
 9028 FORMAT(/,38H TRANSFORMATION OF THE BASIS FUNCTIONS,/)
 9048 FORMAT(7X,15(2X,A4,2X))
 9068 FORMAT(2X,A4,1X,15F8.4)
 9088 FORMAT(1H1)
 9108 FORMAT(/,21X,21HTRANSFORMATION NUMBER,I4,/)
 9128 FORMAT(/)
      END
C*MODULE SCFLIB  *DECK WDISK
      SUBROUTINE WDISK(A1,A2,A3,NDAF,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      DIMENSION A1(2),A2(2),A3(2)
C
      CALL DAWRIT(IDAF,IODA,A1,L2,NDAF  ,0)
      CALL DAWRIT(IDAF,IODA,A2,L2,NDAF+1,0)
      CALL DAWRIT(IDAF,IODA,A3,L2,NDAF+2,0)
C
      RETURN
      END
C*MODULE SCFLIB  *DECK SHLDEND
      SUBROUTINE SHLDEND(SCFTYP,DA,DB,DGVB,DSH,L1,NFO)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL ALPHA,BETA,FMOESP
C
      DIMENSION DA(L1),DB(L1),DGVB(L1,NFO),DSH(*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /FMORUN/ ESPSCF,E0SCF(2),IDAFMO,ICURLAY,ICURFG,JCURFG,
     *                KCURFG,NCURSH,NGAU,ICURPOP,IFMOSTP,MONCOR,NAT1E,
     *                NEEDR,MODRST,NORBPROJ,NUNESP,ISKIPESP,ICURUNT
      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
C
      PARAMETER (ZERO=0.0D+00)
C
      DATA RHF,UHF,ROHF,GVB/8HRHF     ,8HUHF     ,8HROHF    ,8HGVB     /
C
C     --- COMPRESS DIAGONAL DENSITY OVER AOS TO DENSITY OVER SHELLS ---
C        PROCESS ONLY SHELLS FROM THE EXTERNAL MONOMER FOR FMO.
C
      ALPHA = SCFTYP.EQ.RHF  .OR.  SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
      BETA  =                      SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
      FMOESP=NFG.NE.0.AND.NCURSH.NE.0
      ISH0=1
      NEAO=0
      IF(FMOESP) THEN
         ISH0=NCURSH+1
         NEAO=L1
      ENDIF
      DO 240 ISH=ISH0,NSHELL
         MINI = KLOC(ISH)
         MAXI = MINI + KMAX(ISH) - KMIN(ISH)
         DMAX = ZERO
         IF(ALPHA) THEN
            DO I=MINI,MAXI
               DMAX=MAX(DMAX,ABS(DA(I-NEAO)))
            ENDDO
         ENDIF
         IF(BETA) THEN
           DO I=MINI,MAXI
              DMAX=MAX(DMAX,ABS(DB(I-NEAO)))
           ENDDO
         ENDIF
         IF(SCFTYP.EQ.GVB) THEN
            DO 110 IHAM=1,NFO
               DO I=MINI,MAXI
                  DMAX=MAX(DMAX,ABS(DGVB(I-NEAO,IHAM)))
               ENDDO
  110       CONTINUE
         END IF
         DSH(ISH-ISH0+1) = DMAX
  240 CONTINUE
      RETURN
      END
C*MODULE SCFLIB  *DECK PRSELMO
      SUBROUTINE PRSELMO(VEC,EIG,LABMO,NMO,NAO,LDV)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      DIMENSION VEC(LDV,NMO),EIG(NMO),LABMO(NMO)
      PARAMETER (MXATM=2000)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
C     PRINT ORBITALS.
C     NPREO(1)=2,1 IS TREATED AS A SPECIAL CASE (DEFAULT SETTINGS
C     IN $SCF): RATHER THAN PRINTING NOTHING, PRINT EVERYTHING.
C     IF ONE WANTS TO PRINT NOTHING, ANY OTHER "TWISTED" SETTING
C     SHOULD BE USED, E.G. NPREO(1)=3,1.
C
      IF(MASWRK) THEN
      WRITE(IW,9180)
C
      MOST=NPREO(1)
      IF(MOST.LE.0) MOST=MAX(NA+MOST,1)
      MOFI=NPREO(2)
      IF(MOFI.LE.0) MOFI=MIN(NA-MOFI,NMO)
      IF(MOFI.GT.NMO) MOFI=NMO
      NPRMO=MOFI-MOST+1
      IF(NPRMO.GE.1.AND.NPRMO.LE.NMO) THEN
        CALL PREVS(VEC,EIG,LABMO,MOFI,NAO,LDV,MOST)
      ENDIF
C
C     PRINT ENERGIES.
C
      NOST=NPREO(3)
      IF(NOST.LE.0) NOST=MAX(NA+NOST,1)
      NOFI=NPREO(4)
      IF(NOFI.LE.0) NOFI=MIN(NA-NOFI,NMO)
      NOE=NOFI-NOST+1
      IF(NOE.GE.1.AND.NOE.LE.NMO) THEN
        WRITE(IW,9190)
        CALL PRORBEN(EIG,LABMO,NOFI,NOST)
      ENDIF
      ENDIF
      RETURN
 9180 FORMAT(/10X,12(1H-)/10X,12HEIGENVECTORS/10X,12(1H-))
 9190 FORMAT(/10X,11(1H-)/10X,11HEIGENVALUES/10X,11(1H-))
      END
C
C*MODULE MTHLIB  *DECK PRORBEN
      SUBROUTINE PRORBEN(E,LABMO,N,ISTE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION E(N),LABMO(N)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- PRINT OUT A SQUARE MATRIX -----
C     -V- IS -N- ROWS BY -M- COLUMNS, WITH LEADING DIMENSION -NDIM-
C
      IF (MASWRK) THEN
      MAX = 6
      IF (NPRINT .EQ. 6) MAX = 12
      IMAX = ISTE-1
  100 IMIN = IMAX+1
      IMAX = IMAX+MAX
      IF (IMAX .GT. N) IMAX = N
      WRITE (IW,9008)
      WRITE (IW,9028) (I,LABMO(I),I = IMIN,IMAX)
      WRITE (IW,9048) (E(I),I = IMIN,IMAX)
      IF (IMAX .LT. N) GO TO 100
      WRITE (IW,9008)
      END IF
      RETURN
 9008 FORMAT(1X)
 9028 FORMAT(1X,10(2X,I5,1X,A4,1X))
 9048 FORMAT(1X,10F13.6)
      END
