C  4 MAR 08 - HN  - LINMC,MCIN,MCSCF,NEWTON,NTNCAN: ENABLE CISTEP=GMCCI
C  7 DEC 07 - JI  - MCSCF: FIX ORMAS CI VECTOR PRINTING
C 20 AUG 07 - DGF - SYNCHRONISE FMORUN
C 24 MAR 07 - MWS - NEWTON: PRINT ORMAS' CI EXPANSION
C 22 DEC 06 - DGF - SYNCHRONISE FMORUN
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 25 OCT 06 - SPW - MCSCF: INTERFACE TO NEO-MCSCF ADDED
C  7 SEP 06 - DMC - MCSCF: IMPLEMENT FULL SVPE SOLVATION
C 13 MAR 06 - NM  - MCIN,MCSCF,NTNCAN: CHANGES FOR MEX RUN TYPE
C  2 FEB 06 - MWS - MCIN: ALLOW NUMERICAL JACOBI GRADIENT TO RUN
C 17 JAN 06 - MWS - UPDATE ARGS TO TRFMCX
C 14 NOV 05 - DGF,HL - SOLVENT SCF MEMORY CHANGES
C 19 SEP 05 - MWS - CHANGE STRATEGY FOR READING FORS KEYWORD
C 13 JUL 05 - MWS - NTNCAN: MRPT SHOULD FORCE ACTIVE SPACE CANON
C  5 JUL 05 - MWS - ACURCY TEN TIMES TIGHTER IF MPLEVL OR CITYP,
C                   NTNROT,NTNDVD'S CRITERIA DEPEND ON ACURCY
C 27 JUN 05 - MWS - MCSCF: SVP MEMORY ALLOCATION CHANGE
C  1 JUN 05 - MWS - MCSCF: PCM MASWRK I/O, IROOT BUG, MRPT MUST CANON.
C 30 APR 05 - MWS - IMPLEMENT SVP SOLVATION MODEL
C 30 APR 05 - DGF - IMPLEMENT FMO-MCSCF
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL, FRGINF
C  5 FEB 05 - MWS - MCSCF: SET STATE SPECIFIC ENERGY IF REQUESTED
C  7 SEP 04 - MWS - PAD COMMON BLOCK INTFIL
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C  7 MAR 04 - HL  - IEF-PCM/EFP FOR MCSCF
C 16 JAN 04 - GDF - IMPLEMENT DISTRIBUTED MEMORY PARALLEL FULLNR
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  4 NOV 03 - MWS - MCSCF: PASS DDITRF INFORMATION TO ORMAS JOBS
C  3 SEP 03 - MWS - MCIN: ALLOW ACT-ACT OPTIMIZATION IF SOSCF,
C                   MCSCF: CLOSE MOINTS FILE AFTER PARALLEL FULLNR,
C                   NTNDEN: USE JACOBI DIAGONALIZATION
C  3 JUL 03 - MWS - NEWTON: PARALLEL W/O CORE SHOULD USE NTNRGY
C 16 JUN 03 - HL  - MCSCF: PCM GRADIENT CHANGES
C 26 MAR 03 - MWS - MODIFY SQCDF W/O ACURCY CONVERGENCE MESSAGE
C 14 JAN 03 - JI  - ADD HOOKS TO ORMAS PROGRAM
C 12 DEC 02 - MWS - NEWTON: CHANGE DUPLICATED AO LISTS IN ACTIVE FOCK
C  7 AUG 02 - MWS - PRINT ORBITAL PAIR WITH MAX. LAGRANGIAN ASYMMETRY
C  7 AUG 02 - CP,HL - ITERATIVE ISOTROPIC IEF-PCM/EFP
C 20 JUN 02 - MWS - ALLOW USER TO AVOID OUTPUT SUPPRESSION AFTER ITER 1
C 22 MAY 02 - GDF - DISTRIBUTED MEMORY FULLNR ENABLED, DAVIDSON SOLVER
C 17 APR 02 - MWS - MCSCF: CAUTIONARY LABELING OF CI EXPANSION ADDED
C 26 MAR 02 - KRG - USE ABRT CALL INSTEAD OF STOP
C 24 JAN 02 - MWS - MCSCF: WIDEN ITERATION'S ENERGY FIELD
C 16 NOV 01 - MWS - NTNCAN,NTNDEN: CANONICAL ACTIVE ORBS, Q-MATRIX NOS
C  8 OCT 01 - MWS - FIXED ONE OUT OF THREE PCM PARALLELIZATION BUGS
C 19 SEP 01 - MWS - CONVERT MXAOCI PARAMTER TO MXAO
C  6 SEP 01 - HN  - NTNCAN: FIX EIGENVALUE SORTING
C  1 AUG 01 - JI  - VARIOUS CALLS FOR GENERAL CI
C 25 JUN 01 - MWS - ALTER COMMON BLOCK SCFOPT
C 13 JUN 01 - MWS - PARALLEL TEI, NTNCAN: FIX SPHERICAL HARMONICS BUG
C 20 FEB 01 - MWS - NTNCAN: IF FCORE, DON'T CANONICALIZE THE CORE
C 11 OCT 00 - MWS - NTNCAN: FIX LAGRANGIAN ROTATION ORBITAL PHASE BUG
C  6 SEP 00 - BM  - ADDED IEF MODEL
C 15 AUG 00 - MWS - NTNPAR: CHECK RUNS LOOK FOR SYMMETRY LABELS ON DAF
C 28 JUL 00 - MWS - NTNCAN,CANMCV: OBTAIN SYMMETRIC CANONICAL ORBS
C 11 JUN 00 - MWS - MCSCF,NEWTON: ALWAYS ALLOCATE HONDO INTEGRAL BUFFER
C  1 MAY 00 - MWS - NO SYMMETRY OPTION IN TRFOPT COMMON
C 16 FEB 00 - MWS - ADD CONVERGENCE TEST ON SQCDF GOING VERY SMALL
C 21 DEC 99 - MWS - REMOVE C1 DOWNSHIFT FROM TRANSF. BUT NOT AO INTS
C 13 FEB 99 - MWS - NTNONE,NTNRGY: USE GSUMCP TO LIMIT GLOBAL SUM SIZE
C 12 NOV 98 - GDF - CHANGE BIT PACKING TO ISHIFT
C 30 OCT 98 - DGF - NTNROT: FCORE SKIPS CORE ELEMENTS IN LAGRANGIAN
C 26 OCT 98 - MAF - CANMCV; MCSCF; NTNCAN: ALLOW USE OF SPHER. HARMONICS
C 12 SEP 98 - MWS - MCIN: CORRECT CANONC TO BE LOGICAL INPUT
C 13 APR 98 - MWS - MCIN,LINMC,MCSCF: ALLOW USE OF DETERMINANT CI STEP
C 23 MAR 98 - MWS - NTNCAN: SAVE SYMMETRY LABELS ON DAF
C 27 FEB 98 - MWS - NTNCAN: PRINT ALL MCSCF ORBITALS
C  6 JAN 98 - MWS - CHANGE RAOPEN,RACLOS CALLS
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C 21 APR 97 - MWS - NTNPAR: ADDITIONAL MO FREEZING OPTIONS.
C  2 APR 97 - MWS - MCSCF: ZERO DAMPING FACTOR FOR FOCAS AND SOSCF
C 18 MAR 97 - BM  - MCSCF: INCLUDE PCM INTERACTIONS
C 19 FEB 97 - MWS - MCIN: ALLOW FORS TO BE FALSE DURING SOSCF RUNS
C 14 FEB 97 - MWS - MCSCF: ZERO DAMPING AFTER CONVERGENCE OBTAINED
C  7 FEB 97 - MWS - MCSCF: FIX BUG WITH RUNNING EFP IN PARALLEL
C  3 FEB 97 - MWS - MCSCF: FORCE DROPC OFF IF THERE ARE NO CORES
C  8 JAN 97 - GMC - LINMC,MCIN,MCSCF,NEWTON: CHANGES FOR DROPPING MCC-S
C                   NTNRGY,NTNDE1,NTNLAG,NTNHES: NEW FOR FOCK CORE TERMS
C                   NTNONE,NTNRPH: USE SQWRIT AND SQREAD
C 26 NOV 96 - VAG - INCREASE CI AO'S TO 768
C 30 OCT 96 - GMC - MCSCF: AVOID CALL TO NTNDMP EXCEPT FOR FULLNR METHOD
C 23 OCT 96 - SPW - MCSCF: CHANGES FOR SCRF AND IEFP RUNS; ADD MCAODEN
C 29 SEP 96 - GMC - MCSCF: DON'T KILL PARALLEL FOCAS OR SOSCF JOBS
C 10 SEP 96 - MWS - MCIN: KILL JOB IF FCORE IS SELECTED W/O FULLNR
C 13 MAY 96 - VAG - ADD VARIABLE FOR CI TYPE TO SCFOPT COMMON
C 24 MAY 96 - GMC - PROPERLY DOWNSHIFT TO C1 IF NON-ABELIAN SOSCF
C 18 APR 96 - GMC - CHANGES TO CALL FOCAS AND SOSCF CODE
C  5 MAR 96 - MWS - CHANGE SHELL SYMMETRY PACKING COMMONS
C  9 FEB 96 - MWS - MCSCF: FIX FROZEN CORE CONVERGENCE TEST
C 20 NOV 95 - MWS - MCIN: OPTIMIZE ACT-ACT ROTATIONS BY DEFAULT
C 20 OCT 95 - GMC - MODIFY LINSER EXIT CONDITIONS
C  6 OCT 95 - MWS - READ INPUT IN MCIN GROUP
C 28 SEP 95 - GPD - ADD LINEAR SEARCH ROUTINES
C 26 JUL 95 - MWS - NTNCAN: ALLOW USER CONTROL OF ORBITAL PUNCHING
C 12 JUN 95 - MWS - NTNFCK: GLOBAL SUM ADDED FOR DIST. AO INTS
C 21 APR 95 - MWS - USE CORRECT NAME FOR EXCHANGE INTEGRAL SCREENING
C  5 MAR 95 - MWS - MCSCF: CONVERGENCE IN OPT RUNS DEPENDS ON GRADIENT
C 12 NOV 94 - MWS - REMOVE MOST FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C  9 AUG 94 - MWS - NTNFCK: USE DOUBLE LABEL PACKING
C  1 JUN 94 - MWS - MCSCF: CHANGE ARGUMENTS TO TRANSFORMATION CALLS
C 23 APR 94 - MWS - MCSCF: GENERATE EXCHANGE INTEGRALS FOR DIRECT RUNS
C 31 MAR 94 - MWS - ADD A VARIABLE TO END OF MACHSW COMMON
C 28 FEB 94 - MWS - NTNCAN: USE JACOBI DIAGONALIZATION
C 25 JAN 94 - MWS - MCSCF: USE C1 SYM FOR NON-ABELIAN TRANSFORMATIONS
C  4 JAN 94 - MWS - NTNFCK: SYMMETRIZE SKELETON FOCK MATRIX
C 13 DEC 93 - TLW - INCLUDE CHANGES FOR PARALLEL MCSCF
C  9 DEC 93 - MWS - EKTM,NTNDEN,NTNONE: CHANGE DM1'S DAF RECORD NO.
C 16 JUL 93 - MWS - INCREASE MAXIMUM CI ROOTS TO 100
C 24 FEB 93 - KAN,MWS - INCLUDE MORRISON'S EKT OPTION
C  2 APR 92 - MWS,TLW - COMMON ENRGYS MADE PURE FLOATING POINT
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C 25 FEB 92 - MWS - NTNCAN: CHANGE SYMTRZ TO SYMMOS
C  5 FEB 92 - TLW - CANMCV: STOP PRINTING MO LAGRANGIAN
C 11 JAN 92 - TLW - MAKE READ PARALLEL
C 10 JAN 92 - TLW - CHANGE CLOSES TO CALL SEQCLO; CHANGE REWINDS TO
C                   CALL SEQREW; NTNFCK: CALL TO ABRT BEFORE STOP
C 10 JAN 92 - MWS,TLW - BRING OPENFM USING SEQOPN TO THIS MODULE
C  7 JAN 92 - TLW - MAKE WRITES PARALLEL; ADD COMMON PAR
C 15 NOV 91 - MWS - WRITE PRINTOUT DRIVER ROUTINE CANMCV,
C                   DON'T GENERATE CANONICAL ORBITALS DURING GEOM SEARCH
C 31 OCT 91 - JHJ - MCSCF: ADD COMMON /CONV/.
C  7 OCT 91 - MWS - CHANGE UNIX BYTE UNPACKING
C 28 JUN 91 - MWS - TURN OFF DRT PRINTOUT AFTER 1ST GEOMETRY
C 17 APR 91 - MWS - NEWTON: CALL TO ABRT IF THERE IS NOT ENOUGH MEMORY
C 12 SEP 90 - MWS - INTRODUCE MXATM
C  4 JUL 90 - MWS - NTNDEN: USE DRT'S CORE ORBITAL COUNT
C  3 JUN 90 - JAM - CHANGE VAX BYTE PACKING TO FULL 8 BITS
C 10 MAR 90 - MWS - NEWTON: INITIALIZE MEM
C 23 JAN 90 - MWS - NEWTON: FIX BUG IN MEMORY DEALLOC WHEN NO VIRTUALS
C  7 JAN 90 - MWS - ASSIGN ALL MO SYMMETRIES, $DRT REPLACES $DRTINP
C 16 DEC 89 - MWS - DOUBLE NDAR20, AND ALLOW NDAR TO BE READ IN.
C 27 SEP 89 - MWS - CALL FLSHBF AFTER EACH ITERATIONS OUTPUT
C 26 SEP 89 - MWS - USE /CIFILS/ INSTEAD OF DATA TO DEFINE FILES.
C  5 MAY 89 - MWS - ONLY ASSIGN SYMMETRIES TO CORE MCHF ORBITALS.
C  2 MAR 89 - STE - NEWTON: FIX COMPUTATION OF NGOT
C 26 FEB 89 - STE - DELETE FORMAT 9340 FROM NEWTON
C 20 FEB 89 - STE - SPELL LIARC CORRECTLY IN NEWTON
C  1 FEB 89 - MWS - ADD CALL TO SYMTRZ, PREVS IN NTNCAN
C 15 DEC 88 - MWS - CHANGE CALLS TO TRNSF
C 22 OCT 88 - MWS - FIX BUG WITH LEARNING MEMORY NEEDS
C 13 OCT 88 - MWS - IMPLEMENT MICHEL'S HONDO7 MCSCF PROGRAM IN GAMESS
C
C*MODULE MCSCF   *DECK CANMCV
      SUBROUTINE CANMCV(NPFLG8,CANONC,NPUNCH,
     *                  NORB,NCORBS,NVAL,NVIRT,NORBS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXRT=100, MXAO=8192, MXSH=5000,
     *           MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH)
C
      LOGICAL SOME,DBG,PACK2E,CANONC,GOPARR,DSKWRK,MASWRK,
     *        SCHWRZ,DIRTRF
C
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FMCOM / X(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /INT2IC/ NINTIC,ININTIC,NXXIC,LBUFPIC,LIXIC,LABSIX,NINTIX
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      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),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /ORDOPT/ NORDER,NDAR,LDAR,NBOXMX,NWORD,NOMEM,NSQUAR
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
C
      DATA RHF/8HRHF     /
C
C           DRIVER TO GENERATE CANONICAL ORBITALS,
C           AND THEN PRINT NATURAL AND MC-HF ORBITALS.
C
      DBG  = (NPFLG8.GT. 0 .AND. MASWRK)
      SOME = (NPRINT.NE.-5 .AND. MASWRK) .OR. DBG
C
      L0 = NQMT
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
C
      DO 100 I=1,L1
         IA(I) = (I*I-I)/2
  100 CONTINUE
C
      NSH2 = (NSHELL*NSHELL+NSHELL)/2
      CALL BASCHK(LMAX)
                    NANGM =  4
      IF(LMAX.EQ.2) NANGM =  6
      IF(LMAX.EQ.3) NANGM = 10
      IF(LMAX.EQ.4) NANGM = 15
      IF(LMAX.EQ.5) NANGM = 21
      IF(LMAX.EQ.6) NANGM = 28
      MAXG = NANGM**4
C
      MINTMX=NINTMX
      IF(NINTIC.NE.0) MINTMX=0
C     MEMORY IS ALLOCATED ELSEWHERE
C
      CALL VALFM(LOADFM)
      LVEC  = LOADFM+1
      LWRK1 = LVEC  + L3
      LWRK2 = LWRK1 + MAX(L2,L3)
      LFCK  = LWRK2 + MAX(L2,L3)
      LEIG  = LFCK  + L2
      LSCR  = LEIG  + L1
      LIWRK = LSCR  + L1*8
      LOCC  = LIWRK + L1
      LS    = LOCC  + L1
      LQ    = LS    + L2
      LFMO  = LQ    + L3
      LAST  = LFMO  + L2
      IF(DIRTRF) THEN
         LXX    = LAST
         LIXX   = LAST
         LXINTS = LAST
         LGHOND = LXINTS + NSH2
         LDSH   = LGHOND + MAXG
         LDDIJ  = LDSH   + NSH2
         LAST   = LDDIJ  + 49*MXG2
      ELSE
         LXX    = LAST
         LIXX   = LXX   + MINTMX
         LAST   = LIXX  + MINTMX
         LXINTS = LAST
         LGHOND = LAST
         LDSH   = LAST
         LDDIJ  = LAST
      END IF
      IF(NINTIC.NE.0) THEN
         LXX    = LBUFPIC
         LIXX   = LIXIC
      ENDIF
      NEED = LAST-LOADFM-1
      CALL GETFM(NEED)
C
C     ----- CONSTRUCT THE SO-CALLED STANDARD FOCK OPERATOR -----
C
      IF(CANONC) THEN
         LDENS = LWRK1
         IF(DIRTRF) THEN
            IF(SOME) WRITE(IW,9020)
            IDUMMY = 0
            DUMMY  = 0.0D+00
            CALL DAREAD(IDAF,IODA,X(LDENS),L2,16,0)
            CALL VCLR(X(LFCK),1,L2)
            SCHWRZ = ISCHWZ.GT.0
            IF(SCHWRZ) THEN
C                   FOR FMO ORBITAL CONSTRUCTOR EXCHANGE INTEGRALS MAY
C                   NOT HAVE BEEN COMPUTED: DO IT NOW.
C               IF(NORB.EQ.0) CALL EXCHNG(X(LXINTS),X(LGHOND),X(LDDIJ),
C    *                                    NSH2,MAXG,INTTYP)
                CALL DAREAD(IDAF,IODA,X(LXINTS),NSH2,54,0)
                CALL SHLDEN(RHF,X(LDENS),DUMMY,DUMMY,X(LDSH),
     *                      IA,L1,L2,NSH2,1)
            END IF
C              SKELETON FOCK OPERATOR
            CALL TWOEI(RHF,.TRUE.,.FALSE.,.FALSE.,.FALSE.,
     *                 INTTYP,SCHWRZ,NINT,NSCHWZ,L1,L2,
     *                 DUMMY,DUMMY,IDUMMY,1,
     *                 X(LXINTS),NSH2,X(LGHOND),MAXG,X(LDDIJ),
     *                 IA,X(LWRK1),X(LFCK),DUMMY,DUMMY,X(LDSH),
     *                 DUMMY,DUMMY,1)
C              DENSITY SCALING/GLOBAL SUM/SYMMETRIZATION/ADD 1E- TERM
            CALL DSCAL(L2,0.5D+00,X(LFCK),1)
            II=LFCK-1
            DO I=1,L1
               II = II+I
               X(II) = X(II) + X(II)
            ENDDO
            IF(GOPARR) CALL DDI_GSUMF(2020,X(LFCK),L2)
            CALL SYMH(X(LFCK),X(LDENS),IA)
            CALL DAREAD(IDAF,IODA,X(LDENS),L2,11,0)
            CALL VADD(X(LFCK),1,X(LDENS),1,X(LFCK),1,L2)
         ELSE
            IF(SOME) WRITE(IW,9030)
            CALL NTNFCK(X(LDENS),X(LFCK),X(LXX),X(LIXX),
     *                  NINTMX,X(LWRK2),IA,L1,NSQUAR)
         END IF
      END IF
C
      IF(DBG) THEN
         IF(MASWRK) WRITE(IW,*) 'THE STANDARD FOCK OPERATOR IS'
         CALL PRTRI(X(LFCK),L1)
      END IF
C
C     ----- MAKE CANONICAL ORBITALS -----
C     OR, AT LEAST, PRINT OUT THE NO-S AND MO-S
C
      CALL NTNCAN(X(LVEC),X(LWRK2),X(LFMO),X(LFCK),X(LS),X(LQ),
     *            X(LWRK1),X(LOCC),X(LEIG),X(LSCR),X(LIWRK),
     *            NORB,NCORBS,NVAL,NVIRT,NORBS,L0,L1,L2,L3,
     *            ETOT,ENUCR,SOME,DBG,NPUNCH)
C
      CALL RETFM(NEED)
      RETURN
C
 9020 FORMAT(/1X,'FORMING THE "STANDARD FOCK OPERATOR" IN DIRECT AO',
     *           ' INTEGRAL MODE...')
 9030 FORMAT(/1X,'FORMING THE "STANDARD FOCK OPERATOR" USING INTEGRALS',
     *           ' FROM DISK...')
      END
C*MODULE MCSCF   *DECK EKTM
      SUBROUTINE EKTM(NORBS,ACURCY)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DEBUG,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      PARAMETER (ZERO=0.0D+00, HALF=0.5D+00)
C
      DATA DBUG,DBUGME/8HDEBUG   ,8HEKT     /
C
C     ----- ALLOCATE MEMORY FOR EXTENDED KOOPMANS' THEOREM -----
C     THE ORIGINAL VERSION OF THIS ROUTINE WAS WRITTEN BY
C         ROBERT C. MORRISON
C         DEPARTMENT OF CHEMISTRY
C         EAST CAROLINA UNIVERSITY
C         GREENVILLE, NC  27858-4353
C         CHMORRIS @ ECUVM.CIS.ECU.EDU
C     SEE THE PAPERS BY
C         R.C.MORRISON, J.CHEM.PHYS. 1992, 96, 3718-3722
C         R.C.MORRISON, G.LIU, J.COMPUT.CHEM. 1992,13, 1004-1010
C
      DEBUG = (EXETYP.EQ.DBUG  .OR.  EXETYP.EQ.DBUGME)  .AND.  MASWRK
      IF(MASWRK) WRITE(IW,9000)
C
C     NORBS = NUMBER OF OCCUPIED MOLECULAR ORBITALS (CORE + ACT.)
C     NUM   = TOTAL NUMBER OF AOS
C
      M1 =  NORBS
      M2 = (NORBS*NORBS+NORBS)/2
      M3 =  NORBS*NORBS
      L0 =  NQMT
      L1 =  NUM
      L2 = (NUM*NUM+NUM)/2
      L3 =  NUM*NUM
C
      CALL VALFM(LOADFM)
      LEIG  = LOADFM + 1
      LIWRK = LEIG   + M1
      LSCR  = LIWRK  + M1
      LROT  = LSCR   + M1*8
      LDM   = LROT   + M3
      LDMH  = LDM    + M2
      LFOCK = LDMH   + M3
      LVEC  = LFOCK  + M2
      LVEKT = LVEC   + L3
      LS    = LVEKT  + L1*M1
      LAST  = LS     + L2
      NEED  = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
C     READ IN THE FOCK AND DENSITY MATRIX
C  THE FOCK (LAGRANGE MULTIPLIER) MATRIX AND DENSITY MATRIX ARE THE
C  SUM OF ALPHA AND BETA SPIN PARTS.  FOR CLOSED SHELL DIVIDE BY 2.
C  ACTUALLY, THE CODE LETS ANY MULTIPLICTY RUN, UNIMPEDED.
C
      CALL DAREAD(IDAF,IODA,X(LROT),M3,403,0)
      CALL DAREAD(IDAF,IODA,X(LDM) ,M2,68,0)
      CALL DSCAL(M3,HALF,X(LROT),1)
      CALL DSCAL(M2,HALF,X(LDM),1)
C
C     TRANSFER LAGRANGIAN TO LOWER TRIANGULAR STORAGE
C
      IJ=0
      ASYMM=ZERO
      DO 120 I=1,NORBS
         DO 110 J=1,I
            FIJ=X(LROT+(I-1)*NORBS+J-1)
            FJI=X(LROT+(J-1)*NORBS+I-1)
            FTEST = ABS(FIJ-FJI)
            IF(FTEST.GT.ASYMM) ASYMM = FTEST
            IJ=IJ+1
            X(LFOCK+IJ-1)= HALF*(FIJ+FJI)
  110    CONTINUE
  120 CONTINUE
      IF(MASWRK  .AND.  ASYMM.GT.ACURCY) WRITE(IW,9010) ASYMM,ACURCY
C
      CALL EKTORB(X(LEIG),X(LIWRK),X(LSCR),X(LROT),X(LDM),X(LDMH),
     *            X(LFOCK),X(LVEC),X(LVEKT),X(LS),
     *            L0,L1,L2,L3,M1,M2,DEBUG)
      CALL RETFM(NEED)
      RETURN
C
 9000 FORMAT(/1X,'EXTENDED KOOPMANS'' THEOREM I.P. CALCULATION')
 9010 FORMAT(/1X,'*** WARNING *** LAGRANGIAN ASYMMETRY =',1P,E15.6/
     *        1X,' THIS WAVEFUNCTION IS NOT CONVERGED TO',1P,E15.6)
      END
C
C*MODULE MCSCF   *DECK EKTORB
      SUBROUTINE EKTORB(EIG,IWRK,SCR,ROT,DM,DMH,FOCK,VEC,VEKT,S,
     *                  L0,L1,L2,L3,M1,M2,DEBUG)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DEBUG,GOPARR,DSKWRK,MASWRK
C
      DIMENSION EIG(M1),IWRK(M1),SCR(M1,8),ROT(M1,M1),DM(M2),
     *          DMH(M1,M1),FOCK(M2),VEC(L1,M1),VEKT(L1,M1),S(L2)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
C     ----- CALCULATE THE EXTENDED KOOPMANS' THEOREM I.P.'S -----
C     FIRST, THE DENSITY PLAYS THE ROLE OF THE METRIC,
C     SO WE NEED TO COMPUTE ITS RECIPROCAL SQUARE ROOT.
C
      IF(DEBUG) THEN
         WRITE(IW,*) 'DENSITY IN MCSCF MO BASIS IS'
         CALL PRTRI(DM,M1)
         WRITE(IW,*) 'LAGRANGIAN IN MCSCF MO BASIS IS'
         CALL PRTRI(FOCK,M1)
      END IF
C
      CALL GLDIAG(M1,M1,M1,DM,SCR,EIG,ROT,IERR,IWRK)
      IF(IERR.NE.0) THEN
         WRITE(IW,*) 'EKTORB-1: PROBLEM DIAGONALIZING'
         CALL ABRT
      END IF
C
      DO 120 I=1,M1
         EIG(I) = ONE/SQRT(EIG(I))
  120 CONTINUE
C
      DO 160 I=1,M1
         DO 150 J=1,M1
            SUM=ZERO
            DO 140 K=1,M1
               SUM = SUM + ROT(I,K)*EIG(K)*ROT(J,K)
  140       CONTINUE
            DMH(I,J)=SUM
  150    CONTINUE
  160 CONTINUE
C
      IF(DEBUG) THEN
         WRITE(IW,*) 'D**-1/2 IS'
         CALL PRSQ(DMH,M1,M1,M1)
      END IF
C
C     COMPUTE FOCK MATRIX AT -DM- IN THE O.N.B. OF D**-(1/2).
C     DIAGONALIZE IT TO OBTAIN ROTATION TO EKT ORBITALS.
C     BACK TRANSFORM BY THIS ROTATION TO THE MCSCF MO BASIS.
C
      CALL TFTRI(DM,FOCK,DMH,SCR,M1,M1,M1)
      CALL GLDIAG(M1,M1,M1,DM,SCR,EIG,ROT,IERR,IWRK)
      IF(IERR.NE.0) THEN
         IF(MASWRK) WRITE(IW,*) 'EKTORB-2: PROBLEM DIAGONALIZING'
         CALL ABRT
      END IF
      CALL TFSQB(ROT,DMH,SCR,M1,M1,M1)
C
C     OBTAIN AND PRINT THE EKT ORBITALS IN THE AO BASIS
C
      CALL DAREAD(IDAF,IODA,VEC,L3,15,0)
      CALL MRARBR(VEC,L1,L1,M1,ROT,M1,M1,VEKT,L1)
C
C     PRINT EKT "NATURAL TRANSITION ORBITALS" WITH SYMMETRY LABELS
C     NOTE THAT -Q- IS READ INTO -VEC-
C
      CALL DAREAD(IDAF,IODA,VEC,L3,45,0)
      CALL DAREAD(IDAF,IODA,S  ,L2,12,0)
      IF(MASWRK) THEN
         WRITE(IW,9010)
         CALL SYMMOS(IWRK,VEC,S,VEKT,SCR,L0,L1,M1,L1)
         CALL PREVS(VEKT,EIG,IWRK,M1,L1,L1,1)
      END IF
      RETURN
C
 9010 FORMAT(/10X,12(1H-)/10X,'EKT ORBITALS'/10X,12(1H-))
      END
C*MODULE MCSCF   *DECK LINMCM
      SUBROUTINE LINMCM(EMC0,ITER,LCNT,DROPC,NORBS,NCORBS,
     *                  DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *                  DOVVVO,DOVVVV,DOCORE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DROPC,DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *        DOVVVO,DOVVVV,DOCORE
C
      PARAMETER (MXATM=2000)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
C
C        ALLOCATE MEMORY FOR GURU DAS' LINE SEARCH
C
      L1 = NUM
      L3 = L1*L1
C
      CALL VALFM(LOADFM)
      LVEC  = LOADFM + 1
      LVEC0 = LVEC   + L3
      LDVEC = LVEC0  + L3
      LAST  = LDVEC  + L3
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
C        SINCE THE METHOD USES THE PREVIOUS ITERATION'S MO-S,
C        ON THE FIRST ITERATION WE JUST SAVE THEM
C
      CALL DAREAD(IDAF,IODA,X(LVEC),L3,15,0)
      IF(ITER.GE.2 .AND. LCNT.EQ.0) THEN
         CALL DAREAD(IDAF,IODA,X(LVEC0),L3,61,0)
         CALL LINMC(X(LVEC),X(LVEC0),X(LDVEC),L3,EMC0,LCNT,
     *              DROPC,NORBS,NCORBS,
     *              DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *              DOVVVO,DOVVVV,DOCORE)
         CALL DAREAD(IDAF,IODA,X(LVEC) ,L3,15,0)
      END IF
      CALL DAWRIT(IDAF,IODA,X(LVEC),L3,61,0)
      LCNT=LCNT-1
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE MCSCF   *DECK LINMC
      SUBROUTINE LINMC(VEC,VEC0,DVEC,L3,EMC0,LCNT,DROPC,NORBS,NCORBS,
     *                 DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *                 DOVVVO,DOVVVV,DOCORE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,DROPC,SVDSKW,
     *        CANONC,FCORE,FORS,NOCI,EKT,LINSER,
     *        DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,DOVVVO,DOVVVV,DOCORE
      DOUBLE PRECISION METHOD
C
      PARAMETER (MXATM=2000, MXNORO=250)
C
      DIMENSION VEC(L3),VEC0(L3),DVEC(L3)
C
      COMMON /FUNCT / E,EGRAD(3,MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,NORB,
     *                NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      INTEGER         D_OOOO,D_VOOO,D_VVOO,D_VOVO,D_VVVO,D_VVVV,
     *                D_OOOOAB,D_OOOOBB,D_VOOOAB,D_VOOOBA,D_VOOOBB,
     *                D_VVOOAB,D_VVOOBA,D_VVOOBB,D_VOVOAB,D_VOVOBB,
     *                D_U,D_UB,D_E,D_EB
      LOGICAL         NDOOOO,NDVOOO,NDVVOO,NDVOVO,NDVVVO,NDVVVV,NDCORE,
     *                NDVVOOBA,NDVVOOAB,NDVVOOBB,NDVOVOAB,NDVOVOBB,
     *                NDVOOOBA,NDVOOOAB,NDVOOOBB,NDOOOOAB,NDOOOOBB
      COMMON /TRFDMS/ D_OOOO,D_VOOO,D_VVOO,D_VOVO,D_VVVO,D_VVVV,
     *                D_OOOOAB,D_OOOOBB,D_VOOOAB,D_VOOOBA,D_VOOOBB,
     *                D_VVOOAB,D_VVOOBA,D_VVOOBB,D_VOVOAB,D_VOVOBB,
     *                D_U,D_UB,D_E,D_EB,
     *                NDOOOO,NDVOOO,NDVVOO,NDVOVO,NDVVVO,NDVVVV,NDCORE,
     *                NDVVOOBA,NDVVOOAB,NDVVOOBB,NDVOVOAB,NDVOVOBB,
     *                NDVOOOBA,NDVOOOAB,NDVOOOBB,NDOOOOAB,NDOOOOBB
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, HALF=0.5D+00, TEN=1.0D+01)
C
      DATA GUGA,ALDET,GENCI,ORMAS,GMCCI
     *      /8HGUGA    ,8HALDET   ,8HGENCI   ,8HORMAS   ,8HGMCCI   /
C
C     ----- LINE SEARCH FOR ENERGY LOWERING -----
C     THE METHOD IS DUE TO GURU DAS, WHO SEARCHES ONLY FOR
C     A POINT WHERE THE ENERGY GOES DOWN, NOT THE MINIMUM ON
C     THE LINE CONNECTING THE ORBITALS OF THE CURRENT AND
C     PREVIOUS ITERATION.  THE IMPLEMENTATION HERE HAS BEEN
C     ADJUSTED SOMEWHAT FROM THE ORIGINAL.
C
C     ON ENTRY,
C     -E- IS CURRENT ENERGY             (LINE SEARCH STEP 0.0)
C     -EMC0- IS PREVIOUS ITER'S ENERGY  (LINE SEARCH STEP 1.0)
C     -VEC- IS CURRENT ORBITALS
C     -VEC0- IS PREVIOUS ITER'S ORBITALS.
C
      NOPRT=-5
C
      ELIN=ZERO
      LCNT=1
C
C         THE LINE SEARCH IS PERFORMED ONLY IF ENERGY ROSE
C
      IF(E.LE.EMC0) RETURN
      DE = E-EMC0
      IF(ABS(DE).LT.TEN*ENGTOL) RETURN
C
C         THE FIRST TWO POINTS IN THE LINEAR SEARCH ARE ALREADY IN HAND
C
      FACAC=ONE
      IF(MASWRK) THEN
         WRITE(IW,9110)
         WRITE(IW,9120) E,FACAC
         WRITE(IW,9120) EMC0,ZERO
      END IF
      CALL FLSHBF(IW)
C
C         CARRY OUT A SEARCH ON THE LINE CONNECTING THE CURRENT
C         AND PREVIOUS ORBITALS FOR A LOWER ENERGY THAN THAT OF
C         THE PREVIOUS ITERATION.
C
  100 CONTINUE
      IF(E.LT.EMC0) GO TO 700
      DE=E-EMC0
      IF(ABS(DE).LT.ENGTOL) THEN
         LCNT=3
         GO TO 600
      END IF
C
C         GENERATE ORBITALS WITH A FULL BACKWARDS STEP.
C         IT IS LEFT TO -TRFMCX- TO ORTHOGONALIZE
C
      DO 110 I=1,L3
         DVEC(I)=VEC(I)-VEC0(I)
  110 CONTINUE
      DO 120 I=1,L3
         VEC(I)=VEC0(I)-DVEC(I)
  120 CONTINUE
      FACAC=-FACAC
C
C         OBTAIN THEIR ENERGY
C
C        NORMALLY THE DDI MATRICES ARE DESTROYED BY CALLING -NEWTON-.
C        THE DISTRIBUTED MATRICES OF PREVIOUS ORBS MUST BE EXPLICITLY
C        DESTROYED NOW, SINCE NO FULLNR ORBITAL UPDATE HAS OCCURED.
C        (DESTROYS IN REVERSE ORDER TO CREATES, SEE TRANDDI)
C
      IF(DDITRF) THEN
         CALL DDI_DESTROY(D_VVOO)
         CALL DDI_DESTROY(D_OOOO)
         CALL DDI_DESTROY(D_VOOO)
         CALL DDI_DESTROY(D_VOVO)
      END IF
C
      CALL DAWRIT(IDAF,IODA,VEC,L3,15,0)
      SVDSKW=DSKWRK
      DSKWRK=.FALSE.
      IF (DROPC) THEN
         CALL TRFMCX(NOPRT,NCORBS,NORBS,NORB,.TRUE.,.FALSE.,
     *               DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *               DOVVVO,DOVVVV,DOCORE)
      ELSE
         CALL TRFMCX(NOPRT,0,0,0,.FALSE.,.FALSE.,
     *               DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *               DOVVVO,DOVVVV,DOCORE)
      END IF
      DSKWRK=SVDSKW
C
      IF(CISTEP.EQ.GUGA) THEN
         SVDSKW=DSKWRK
         DSKWRK=.TRUE.
         CALL GUGSRT(NOPRT,DDITRF)
         CALL GUGAEM(NOPRT)
         CALL GUGADG(NOPRT)
         CALL GUG2DM(NOPRT)
         DSKWRK=SVDSKW
      END IF
      IF(CISTEP.EQ.ALDET) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL DETFCI(NOPRT,.TRUE.,DDITRF)
         CALL DETDM2(NOPRT)
         DSKWRK=SVDSKW
      END IF
      IF(CISTEP.EQ.GENCI) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL DETGCI(NOPRT,.TRUE.,DDITRF)
         CALL GCIDM2(NOPRT)
         DSKWRK=SVDSKW
      ENDIF
      IF(CISTEP.EQ.ORMAS) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL DEFCCI(NOPRT,.TRUE.,DDITRF,0,1,0,-5)
         DSKWRK=SVDSKW
      END IF
      IF(CISTEP.EQ.GMCCI) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL OS_GMCCI_SOLVE(DDITRF)
         DSKWRK=SVDSKW
      END IF
C
      IF(MASWRK) WRITE(IW,9120) E,FACAC
      CALL FLSHBF(IW)
C
      DE1=E-EMC0
      IF(DE1.LT.ZERO) GO TO 700
      IF(ABS(DE1).LT.ENGTOL) THEN
         LCNT=3
         GO TO 600
      END IF
C
C        DO A PARABOLIC FIT, GENERATE ORBITALS AT ITS MINIMUM
C
      FAC=HALF*(DE1-DE)/(DE1+DE)
      FACAC=FACAC*FAC
      DO 210 I=1,L3
         VEC(I)=VEC0(I)+FAC*DVEC(I)
  210 CONTINUE
C
C        EVALUATE THE ENERGY OF THIS ORBITAL SET
C
      CALL DAWRIT(IDAF,IODA,VEC,L3,15,0)
C
      IF(DDITRF) THEN
         CALL DDI_DESTROY(D_VVOO)
         CALL DDI_DESTROY(D_OOOO)
         CALL DDI_DESTROY(D_VOOO)
         CALL DDI_DESTROY(D_VOVO)
      END IF
C
      SVDSKW=DSKWRK
      DSKWRK=.FALSE.
      IF (DROPC) THEN
         CALL TRFMCX(NOPRT,NCORBS,NORBS,NORB,.TRUE.,.FALSE.,
     *               DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *               DOVVVO,DOVVVV,DOCORE)
      ELSE
         CALL TRFMCX(NOPRT,0,0,0,.FALSE.,.FALSE.,
     *               DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *               DOVVVO,DOVVVV,DOCORE)
      END IF
      DSKWRK=SVDSKW
C
      IF(CISTEP.EQ.GUGA) THEN
         SVDSKW=DSKWRK
         DSKWRK=.TRUE.
         CALL GUGSRT(NOPRT,DDITRF)
         CALL GUGAEM(NOPRT)
         CALL GUGADG(NOPRT)
         CALL GUG2DM(NOPRT)
         DSKWRK=SVDSKW
      END IF
C
      IF(CISTEP.EQ.ALDET) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL DETFCI(NOPRT,.TRUE.,DDITRF)
         CALL DETDM2(NOPRT)
         DSKWRK=SVDSKW
      END IF
C
      IF (CISTEP.EQ.GENCI) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL DETGCI(NOPRT,.TRUE.,DDITRF)
         CALL GCIDM2(NOPRT)
         DSKWRK=SVDSKW
      END IF
C
      IF(CISTEP.EQ.ORMAS) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL DEFCCI(NOPRT,.TRUE.,DDITRF,0,1,0,-5)
         DSKWRK=SVDSKW
      END IF
C
      IF(MASWRK) WRITE(IW,9120) E,FACAC
      CALL FLSHBF(IW)
      DLIN=E-ELIN
      IF(ABS(DLIN).LT.ENGTOL) THEN
         LCNT=3
         GO TO 800
      END IF
      IF(CISTEP.EQ.GMCCI) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL OS_GMCCI_SOLVE(DDITRF)
         DSKWRK=SVDSKW
      END IF
C
C        NO LOWERING FOUND, CONTINUE THE SEARCHING
C
      ELIN=E
      GO TO 100
C
C        ENERGY IS NOT CHANGING SIGNIFICANTLY, THEREFORE EXIT.
C        IT IS NOT LIKELY THAT WE HIT THIS EXIT POINT, I HOPE!
C        WE LET THE NEWTON-RAPHSON CONTINUE WITH THE ORBITALS
C        AND DENSITY JUST GENERATED, MAYBE THE N-R WILL SEEM
C        TO BE CONVERGED AS WELL.
C        DO NOT USE LINSER NEXT TWO ITERATIONS
C
  600 CONTINUE
      IF(MASWRK) WRITE(IW,9130) ENGTOL,FACAC
      CALL FLSHBF(IW)
      RETURN
C
C        WE'VE FOUND ORBITALS WITH ENERGY LOWER THAN PREVIOUS ITER,
C        SO WE'LL RETURN, AND LET THEM GO INTO A NEWTON-RAPHSON STEP.
C        WE'RE READY TO PROCEED DIRECTLY INTO THAT STEP, AS THE 2ND
C        ORDER DENSITY HAS BEEN MADE, AND IT ALLOWS A FULL VARIATION
C        OF THE ORBITALS TO OCCUR.
C
  700 CONTINUE
      IF(MASWRK) WRITE(IW,9140)
      CALL FLSHBF(IW)
      RETURN
C
C      ENERGY DOESN'T CHANGE BETWEEN TWO SUCCESSIVE PARABOLA SEARCHES,
C      RETURN TO NEWTON-RAPHSON AND PERFORM TWO ITERATIONS
C      WITHOUT CALLING LINMC ROUTINE
C
  800 CONTINUE
      IF(MASWRK) WRITE(IW,9150)
      CALL FLSHBF(IW)
      RETURN
C
 9110 FORMAT(1X,'ENERGY HAS RISEN FROM PREVIOUS ITERATION,',
     *          ' INITIATING LINE SEARCH...')
 9120 FORMAT(1X,'E=',F19.9,' AT LINE SEARCH COORD=',F8.3)
 9130 FORMAT(1X,'LINE SEARCH HAS CONVERGED TO',1P,E10.1,0P,
     *          ' AT FAC=',F8.3,/
     *          ' DO NOT TURN ON LINE SEARCH NEXT TWO ITERATIONS')
 9140 FORMAT(1X,'LINE SEARCH HAS RESULTED IN ENERGY LOWERING,',
     *          ' RETURNING TO NEWTON-RAPHSON...')
 9150 FORMAT(1X,'ENERGY DOES NOT CHANGE WITHIN LINE SEARCH,',
     *          ' RETURN TO NEWTON-RAPHSON',/' AND PERFORM TWO',
     *          ' ITERATIONS WITHOUT CALLING LINE SEARCH ROUTINE')
      END
C*MODULE MCSCF   *DECK MCAODEN
      SUBROUTINE MCAODEN(V,GMN,GIJ,NUM,NORBS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      DIMENSION GIJ(*),GMN(*),V(NUM,NUM)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- FORM TOTAL DENSITY MATRIX IN -AO- BASIS -----
C
      N3=NUM*NUM
      L2=(NORBS*NORBS+NORBS)/2
      CALL DAREAD(IDAF,IODA,GIJ,L2,68,0)
      CALL DAREAD(IDAF,IODA,V,N3,15,0)
      MN=0
      DO 740 M=1,NUM
         DO 730 N=1,M
            DUM=ZERO
            DO 720 I=1,NORBS
               DO 710 J=1,NORBS
                  II=MAX(I,J)
                  JJ=MIN(I,J)
                  IJ=IA(II)+JJ
                  DNORB=GIJ(IJ)
                  DUM=DUM+DNORB*V(M,I)*V(N,J)
  710          CONTINUE
  720       CONTINUE
            MN=MN+1
            GMN(MN)=DUM
  730    CONTINUE
  740 CONTINUE
      N2=(NUM*(NUM+1))/2
      CALL DAWRIT(IDAF,IODA,GMN,N2,16,0)
      RETURN
      END
C*MODULE MCSCF   *DECK MCIN
      SUBROUTINE MCIN
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION METHOD
C
      LOGICAL CANONC,FCORE,FORS,NOCI,EKT,LINSER,OK,
     *        GOPARR,DSKWRK,MASWRK,FOCAS,SOSCF,FULLNR,QUD,DROPC,
     *        JACOBI,DIRTRF,LFORS,FDIRCT,QCORR
C
      PARAMETER (MXNORO=250, MXATM=2000)
C
      COMMON /JACOBI/ JACOBI,NJAOR,ELAST,ISTAT
      COMMON /CASOPT/ CASHFT,CASDII,NRMCAS,FOCAS,SOSCF,DROPC
      COMMON /FCCWFN/ NSPACE,MSTA(51),MNUM(51),MINI(51),MAXI(51),
     *                IAMI(51),IAMA(51),IBMI(51),IBMA(51),IDIM(51),
     *                LBST(51),NREF0,FDIRCT,QCORR,C0SQ
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INT2IC/ NINTIC,ININTIC,NXXIC,LBUFPIC,LIXIC,LABSIX,NINTIX
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,NORB,
     *                NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /QUDMC / QUDTHR,QUD
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00)
C
C        SET UP NAMELIST SIMULATION FOR $MCSCF GROUP
C
      PARAMETER (NNAM=31)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA RNAM/8HMCSCF   /
      DATA QNAM/8HMETHOD  ,8HACURCY  ,8HENGTOL  ,8HDAMP    ,
     *          8HCANONC  ,8HFCORE   ,8HFORS    ,
     *          8HNOCI    ,8HCISTEP  ,8HMAXIT   ,8HMICIT   ,
     *          8HNWORD   ,8HNORB    ,8HNOROT   ,
     *          8HEKT     ,8HNPUNCH  ,8HLINSER  ,
     *          8HFOCAS   ,8HCASHFT  ,8HCASDII  ,8HNRMCAS  ,
     *          8HSOSCF   ,8HFULLNR  ,8HNPFLG   ,8HNOFO    ,
     *          8HQUAD    ,8HQUDTHR  ,8HDROPC   ,8HMOFRZ   ,
     *          8HJACOBI  ,8HMCFMO   /
      DATA KQNAM/5,3,3,3,  0,0,0,    0,5,1,1,   1,1,-1,   0,1,0,
     *           0,3,3,1,  0,0,101,1,   0,3,0,151,   0,1/
C
      DATA DM2/8HDM2     /
      DATA ENGWRD/8HENERGY  /
      DATA GUGA,ALDET,GENCI,ORMAS,GMCCI
     *    /8HGUGA    ,8HALDET   ,8HGENCI   ,8HORMAS   ,8HGMCCI   /
      DATA RNONE/8HNONE    /
C
      EQUIVALENCE (LFORS,IFORS)
C
C     ----- READ NAMELIST $MCSCF, STORE IN /MCIN/ COMMON  -----
C
      METHOD=DM2
      ACURCY=1.0D-05
      IF(MPLEVL.GT.0)    ACURCY=1.0D-06
      IF(CITYP.NE.RNONE) ACURCY=1.0D-06
      IF(NFG.NE.0)       ACURCY=1.0D-06
      ENGTOL=1.0D+17
      DAMP=ZERO
C     THE ORIGIN OF THE PROBLEMS IN CANONICALISING THE ORBITALS IN
C     FMO-MCSCF IS NOT CLEAR. IT SEEMS LIKELY THAT IT IS NOT A BUG.
C     IT MANY CASES CANONICALISATION WORKS THOUGH.
      CANONC=.TRUE.
      IF(NFG.NE.0) CANONC=.FALSE.
      FCORE=.FALSE.
      IFORS= -278   ! MEANT TO BE TOTALLY STRANGE, EQUIVALENT TO -LFORS-
      NOCI=.FALSE.
      EKT=.FALSE.
      NPUNCH=2
      LINSER=.FALSE.
C
      FOCAS =.FALSE.
      SOSCF =.FALSE.
      FULLNR=.FALSE.
      QUD   =.FALSE.
      JACOBI=.FALSE.
C
      CISTEP=ALDET
C
      MAXIT=0
      MICIT=0
      CASHFT=1.0D+00
      CASDII=5.0D-02
      NRMCAS=1
      QUDTHR=5.0D-02
      DROPC=.TRUE.
      NWORD=0
      NORB=NUM
      KQNAM(14) = 1 + 2*MXNORO*10
      DO 10 I=1,MXNORO
         NOROT(1,I)=0
         NOROT(2,I)=0
   10 CONTINUE
      DO 11 I=1,15
         MOFRZ(I) = 0
   11 CONTINUE
      DO 13 I=1,10
         NPFLG(I) = 0
   13 CONTINUE
      NOFO=1
      MCFMO=0
      IF(NFG.NE.0) MCFMO=1
C
      JRET = 0
      CALL NAMEIO(IR,JRET,RNAM,NNAM,QNAM,KQNAM,
     *            METHOD,ACURCY,ENGTOL,DAMP,CANONC,FCORE,
     *            LFORS,NOCI,CISTEP,MAXIT,MICIT,NWORD,NORB,NOROT,
     *            EKT,NPUNCH,LINSER,FOCAS,CASHFT,CASDII,NRMCAS,
     *            SOSCF,FULLNR,NPFLG,NOFO,QUD,QUDTHR,DROPC,MOFRZ,
     *            JACOBI,MCFMO,
     *            0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,    0,0,0,0,0,   0,0,0,0,0)
      IF(JRET.EQ.2) CALL ABRT
C
C        PICK THE STRONGEST CONVERGER SELECTED, SWITCH OTHERS OFF
C
      IF(QUD) THEN
         FOCAS =.FALSE.
         SOSCF =.FALSE.
         FULLNR=.FALSE.
         JACOBI=.FALSE.
      END IF
      IF(JACOBI) THEN
         FOCAS =.FALSE.
         SOSCF =.FALSE.
         FULLNR=.FALSE.
      END IF
      IF(FULLNR) THEN
         FOCAS =.FALSE.
         SOSCF =.FALSE.
      END IF
      IF(SOSCF) THEN
         FOCAS =.FALSE.
      END IF
C
C        IF NO CONVERGER WAS SELECTED, PICK THE DEFAULT
C
      IF( .NOT.FOCAS   .AND.  .NOT.SOSCF .AND.
     *    .NOT.FULLNR  .AND.  .NOT.QUD   .AND.
     *    .NOT.JACOBI) SOSCF=.TRUE.
C
      IF(ENGTOL.EQ.1.0D+17) THEN
         IF(JACOBI) THEN
            ENGTOL=1.0D-08
         ELSE
            ENGTOL=1.0D-10
         END IF
      END IF
C
      IF (FOCAS  .AND. MICIT.EQ.0) MICIT=5
      IF (SOSCF  .AND. MICIT.EQ.0) MICIT=5
      IF (FULLNR .AND. MICIT.EQ.0) MICIT=1
      IF (QUD    .AND. MICIT.EQ.0) MICIT=1
      IF (JACOBI .AND. MICIT.EQ.0) MICIT=100
      IF (FOCAS  .AND. MAXIT.EQ.0) MAXIT=100
      IF (SOSCF  .AND. MAXIT.EQ.0) MAXIT=60
      IF (FULLNR .AND. MAXIT.EQ.0) MAXIT=30
      IF (QUD    .AND. MAXIT.EQ.0) MAXIT=30
      IF (JACOBI .AND. MAXIT.EQ.0) MAXIT=100
C
C        FORCE CERTAIN OPTIONS DEPENDING ON MAIN OPTIONS
C
C        THE -FULLNR- AND -QUAD- CONVERGERS HAVE ROBUST ACTIVE-ACTIVE
C        ROTATION OPTIMIZATION, BUT -SOSCF- HAS LESS ROBUST CODE FOR
C        THESE.  -FOCAS- LACKS ANY PROVISION FOR A-A AT ALL!  IF WE
C        USE -SOSCF- ON A COMPLETE ACTIVE SPACE, WE PREFER TO SKIP
C        EXPLICITLY OPTIMIZING ACTIVE-ACTIVE UNLESS THE USER REQUESTS
C        DOING SO.
C        NOTE: -LFORS- AND -IFORS- ARE EQUIVALENT STORAGE, AND SO THE
C        VALUE -IFORS- IS A TRICK TO SEE IF THE VALUE WAS SKIPPED IN
C        THE USER'S INPUT FILE.
C
      FORS=.FALSE.
      IF(SOSCF  .AND.  IFORS.EQ.-278) THEN
         IF(CISTEP.EQ.GUGA  .OR.  CISTEP.EQ.ALDET) THEN
            FORS=.TRUE.
         ELSE
            FORS=.FALSE.
         END IF
      END IF
C
C         IF THE USER INPUT SOMETHING, BELIEVE THEM (EXCEPT FOCAS)
      IF(IFORS.NE.-278) FORS=LFORS
      IF(FOCAS)         FORS=.TRUE.
C
      IF(FORS) THEN
         IF((CISTEP.EQ.ORMAS.AND.NSPACE.GT.1) .OR.
     *       CISTEP.EQ.GENCI                  .OR. 
     *       CISTEP.EQ.GMCCI) THEN
            IF(MASWRK) WRITE(IW,9330) CISTEP
            CALL ABRT
         END IF
      END IF
C
      IF(FOCAS .OR. SOSCF .OR. QUD .OR. JACOBI) LINSER=.FALSE.
      IF(FOCAS .OR. SOSCF .OR.          JACOBI)   DROPC=.TRUE.
      IF((FULLNR .OR. QUD) .AND. METHOD .NE. DM2) DROPC=.FALSE.
C
      IF(FCORE) THEN
         CANONC=.FALSE.
         IF(MPLEVL.NE.0) THEN
            IF(MASWRK) WRITE(IW,9310)
            CALL ABRT
         END IF
      END IF
C
C        FORCE A DIRECT TRANSFORMATION IF RUNNING IN PARALLEL
C
      IF (GOPARR  .AND.  FULLNR) THEN
         IF(.NOT.DIRTRF  .AND.  MASWRK) WRITE(IW,9210)
         DIRTRF=.TRUE.
      END IF
C
C        ECHO ALL INPUT EXCEPT NPFLG/NOCI
C
      IF(MASWRK) THEN
         WRITE(IW,9100)
         WRITE(IW,9110) FOCAS,SOSCF,FULLNR,QUD,JACOBI
         WRITE(IW,9115) CISTEP
         WRITE(IW,9120) MAXIT,MICIT,ACURCY,DAMP,CANONC,ENGTOL,
     *                  EKT,NPUNCH,NWORD
         IF(FOCAS)  WRITE(IW,9130) CASHFT,CASDII,NRMCAS
         IF(SOSCF)  WRITE(IW,9140) FORS,NOFO
         IF(FULLNR  .OR.  QUD) THEN
            WRITE(IW,9150) METHOD,FORS,DROPC,LINSER,FCORE,NORB
         END IF
         IF(JACOBI) WRITE(IW,9160) FORS,NORB
         IF(QUD)    WRITE(IW,9170) QUDTHR
      END IF
C
      NNORO = 0
      DO 200 I=1,MXNORO
         IF(NOROT(1,I).NE.0) NNORO=NNORO+1
  200 CONTINUE
      IF(NNORO.GT.0) THEN
         IF(MASWRK) WRITE(IW,9180)
         DO 210 I=1,NNORO
            IF(MASWRK) WRITE(IW,'(10(1X,2I4))') NOROT(1,I),NOROT(2,I)
  210    CONTINUE
      END IF
C
      NFRZ = 0
      DO 230 I=1,15
         IF(MOFRZ(I).NE.0) NFRZ=NFRZ+1
  230 CONTINUE
      IF(NFRZ.GT.0) THEN
         IF(MASWRK) WRITE(IW,9190) (MOFRZ(I),I=1,NFRZ)
      END IF
      IF(NFG.NE.0) THEN
         IF(MASWRK) WRITE(IW,9200) MCFMO
      END IF
C
C         ERROR CHECKING...
C
      NERR = 0
C
      IF(ENGTOL.LT.0.0D+00  .OR.  ENGTOL.GT.1.0D-01) THEN
         IF(MASWRK) WRITE (IW,*) 'ENGTOL IS CRAZY,',ENGTOL
         NERR=NERR+1
      END IF
C
      OK=.FALSE.
      IF(CISTEP.EQ.GUGA)  OK=.TRUE.
      IF(CISTEP.EQ.ALDET) OK=.TRUE.
      IF(CISTEP.EQ.GENCI) OK=.TRUE.
      IF(CISTEP.EQ.ORMAS) OK=.TRUE.
      IF(CISTEP.EQ.GMCCI) OK=.TRUE.
      IF(.NOT.OK) THEN
         IF(MASWRK) WRITE(IW,9220) CISTEP
         NERR=NERR+1
      END IF
C
      IF((NORB.NE.NUM  .OR.  FCORE  .OR.  NFRZ.GT.0  .OR.  NNORO.GT.0)
     *          .AND. RUNTYP.NE.ENGWRD) THEN
         IF (MASWRK) WRITE(IW,9230)
         NERR = NERR+1
      END IF
C
      IF((FOCAS .OR. SOSCF) .AND.
     *   (FCORE .OR. NNORO.GT.0 .OR. NFRZ.GT.0)) THEN
         IF(MASWRK) WRITE(IW,9240)
         NERR=NERR+1
      END IF
C          THE DETERMINANT CI CODES -MUST- HAVE CORE INTEGRALS REMOVED.
      IF(.NOT.DROPC  .AND.  (CISTEP.EQ.ALDET .OR.
     *                       CISTEP.EQ.GENCI .OR.
     *                       CISTEP.EQ.ORMAS .OR.
     *                       CISTEP.EQ.GMCCI))   THEN
         IF(MASWRK) WRITE(IW,9250)
         NERR=NERR+1
      END IF
C
C        THE JACOBI CONVERGER DOES NOT GENERATE THE ORBITAL LAGRANGIAN
C        THE JACOBI CONVERGER IS NOT ENABLED FOR PARALLEL EXECUTION
C
      IF(JACOBI) THEN
         CALL DERCHK(NDER)
         IF(NDER.GT.0  .AND.  NGLEVL.EQ.0) THEN
            IF(MASWRK) WRITE(IW,9260) RUNTYP
            NERR=NERR+1
         END IF
         IF(GOPARR) THEN
            IF(MASWRK) WRITE(IW,9270)
            NERR=NERR+1
         END IF
      END IF
C
C        THE QUAD CONVERGER RELIES ON THE GUGA CI PACKAGE.
C        THE QUAD CONVERGER IS NOT ENABLED FOR PARALLEL EXECUTION,
C        SEE NOTES IN -MCQUD- ABOUT WHAT NEEDS TO BE DONE TO FINISH.
C
      IF(QUD) THEN
         IF(CISTEP.NE.GUGA) THEN
            IF(MASWRK) WRITE(IW,9280)
C
            NERR=NERR+1
         END IF
         IF(GOPARR) THEN
            IF(MASWRK) WRITE(IW,9290)
            NERR=NERR+1
         END IF
      END IF
C     IN-CORE MCSCF RUNS ONLY SEQUENTIALLY.
      IF(NINTIC.NE.0.AND.GOPARR) THEN
         IF(MASWRK) WRITE(IW,9300)
         NERR=NERR+1
      END IF
C
      IF(NERR.GT.0) THEN
         IF(MASWRK) WRITE(IW,9400)
         CALL ABRT
      END IF
      RETURN
C
 9100 FORMAT(/10X,22(1H-)/10X,'MCSCF INPUT PARAMETERS'/10X,22(1H-))
 9110 FORMAT(1X,'CONVERGER SELECTION: FOCAS = ',L1,2X,'SOSCF = ',L1,
     *       2X,'FULLNR = ',L1,2X,'QUD = ',L1,2X,'JACOBI = ',L1)
 9115 FORMAT(1X,'SECULAR EQUATION METHOD CISTEP = ',A8)
 9120 FORMAT(1X,'--- GENERAL INPUT OPTIONS:'/
     *       1X,'MAXIT =',I8,2X,'MICIT =',I8,2X,'ACURCY=',1P,E10.3/
     *       1X,'DAMP  =',F8.3,2X,'CANONC=',L8,2X,'ENGTOL=',1P,E10.3/
     *       1X,'EKT   =',L8,2X,'NPUNCH=',I8,2X,'NWORD =',I10)
 9130 FORMAT(1X,'--- INPUT FOR FOCAS CONVERGER:'/
     *       1X,'CASHFT=',1P,E10.3,4X,'CASDII=',1P,E10.3,
     *       4X,'NRMCAS=',I8)
 9140 FORMAT(1X,'--- INPUT FOR SOSCF CONVERGER:'/
     *       1X,'FORS  =',L8,2X,'NOFO  =',I8)
 9150 FORMAT(1X,'--- INPUT FOR FULLNR CONVERGER:'/
     *       1X,'METHOD=',A8,2X,'FORS  =',L8,2X,'DROPC =',L8/
     *       1X,'LINSER=',L8,2X,'FCORE =',L8,2X,'NORB  =',I8)
 9160 FORMAT(1X,'--- INPUT FOR JACOBI CONVERGER:'/
     *       1X,'FORS  =',L8,2X,'NORB  =',I8)
 9170 FORMAT(1X,'--- INPUT FOR QUAD CONVERGER:'/
     *       1X,'QUDTHR=',1P,E10.3)
 9180 FORMAT(1X,'PAIRS OF FROZEN ORBITAL ROTATIONS ARE')
 9190 FORMAT(1X,'FROZEN ORBITALS ARE:',15I4)
 9200 FORMAT(1X,'--- INPUT FOR FMO-MCSCF:'/1X,'MCFMO =',I8)
C
 9210 FORMAT(1X,'INTEGRAL TRANSFORMATION OPTION -DIRTRF- FORCED TRUE,'/
     *       1X,'DUE TO THE CHOICE OF -FULLNR- IN A PARALLEL JOB.')
 9220 FORMAT(1X,'*** ERROR, CISTEP=',A8,' IS UNRECOGNIZED')
 9230 FORMAT(1X,'*** ERROR, RUNTYP MUST BE ENERGY IF FCORE CHOSEN,',
     *          ' OR NORB.NE.NUM,'/
     *       1X,'OR WHEN THE -MOFRZ- OR -NOROT- OPTION IS SELECTED,'/
     *       1X,'AS THESE OPTIONS PRECLUDE GRADIENT COMPUTATION.')
 9240 FORMAT(1X,'*** ERROR, -FCORE,MOFRZ,NOROT- MUST BE USED WITH'/
     *          ' ONLY THE -FULLNR- OR -JACOBI- CONVERGERS.')
 9250 FORMAT(1X,'*** ERROR, ONLY CISTEP=GUGA ALLOWS DROPC=.FALSE.'/
     *       1X,'PLEASE CHANGE EITHER DROPC OR CISTEP, AND RESUBMIT.')
 9260 FORMAT(1X,'*** ERROR, JACOBI MCSCF OPTIMIZATION HAS NO GRADIENT.'/
     *       1X,'RUNTYP=',A8,' IS IMPOSSIBLE WITH JACOBI CONVERGER.')
 9270 FORMAT(1X,'*** ERROR, THE JACOBI CONVERGER DOES NOT RUN IN',
     *          ' PARALLEL.')
 9280 FORMAT(1X,'*** ERROR, THE QUAD CONVERGER REQUIRES CISTEP=GUGA',
     *          ' AND A $DRT')
 9290 FORMAT(1X,'*** ERROR, THE QUAD CONVERGER DOES NOT RUN IN',
     *          ' PARALLEL.')
 9300 FORMAT(1X,'*** ERROR, IN-CORE MCSCF DOES NOT RUN IN PARALLEL.')
 9310 FORMAT(1X,'FROZEN ORBITALS AND MULTIREFERENCE PT REQUIRE',
     *          ' SEPARATE RUNS.')
 9330 FORMAT(/1X,'*** ERROR ***, CISTEP=',A8/
     *        1X,'FORS=.TRUE. REQUIRES CHOOSING A FULL ACTIVE SPACE,'/
     *        1X,'RATHER THAN CISTEP=ORMAS, GENCI, GMCCI')
 9400 FORMAT(/1X,'ONE OR MORE MISTAKES WERE FOUND IN YOUR $MCSCF',
     *          ' INPUT.'/
     *       1X,'PLEASE FIX THEM, AND TRY YOUR RUN AGAIN.')
      END
C*MODULE MCSCF   *DECK MCSCF
      SUBROUTINE MCSCF
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION METHOD,DETNEO
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH,
     *           MXATM=2000, MXRT=100, MXFRG=1050)
C
      LOGICAL MICCVG,CVGED,CVGING,DMPING,FINAL,MCNOCI,
     *        SOME,GOPARR,DSKWRK,MASWRK,SVDSKW,DIRTRF,PACK2E,
     *        CVGSVP,CVGSVPNR,CVGSVPDN,CVGSVPDN1
      LOGICAL CANONC,FCORE,FORS,NOCI,EKT,LINSER
      LOGICAL INITAL,QUD,FULLNR,SOSCF,FOCAS,IZRFLS,MCDENAO,MINMEM
      LOGICAL DROPC,DOFOCK,DOEXCH
      LOGICAL DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,DOVVVO,DOVVVV,DOCORE
      LOGICAL JACOBI,ECON,NOTRF
      LOGICAL IPCFP
C
      PARAMETER (MXNORO=250, MDIIS=19, MXAO=8192, MXCAV=5810,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /CASOPT/ CASHFT,CASDII,NRMCAS,FOCAS,SOSCF,DROPC
      COMMON /CIFILS/ NFT11,NFT12,NFT13,NFT14,NFT15,NFT16,IDAF20,NEMMX
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,SDET,SZDET,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORBDT,
     *                NADET,NBDET,KDET,KSTDET,IROOT,IPURES,MAXW1,NITDET,
     *                MAXP,NCIDET,IGPDET,KSTSYM,NFTGCI
      COMMON /DIICAS/ NCALL,NDIIS,IDARE(MDIIS*2),IDARF(MDIIS*2)
      COMMON /DIIFIL/ IDAF30,NAV30,IODA30(4*MDIIS+5)
      COMMON /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FMCOM / X(1)
      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 /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FUNCT / E,EGRAD(3,MXATM)
      COMMON /GUGWFN/ NFZC,NMCC,NDOC,NAOS,NBOS,NALP,NVAL,NEXT,NFZV,
     *                IFORS,IEXCIT,ICICI,NOIRR
      COMMON /INDDIP/ TMUX,TMUY,TMUZ,TMUXD,TMUYD,TMUZD,MINMEM
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /JACOBI/ JACOBI,NJAOR,ELAST,ISTAT
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,NORB,
     *                NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /MEXOPT/ MEXSKPGES,MEXSTATE
      COMMON /NEOJOB/ NEORUN,NELERM
      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 /NTNOPT/ ENERGY,ENERG0,DEMAX,SQCDF,ITER,MICRO,NOTRF
      COMMON /ORBSET/ NORBMX,NORBS,NCORBS,NLEVS,NNA,NNB,NNC,NSYM,MSYM,
     *                IDOCC,IVAL,IMCC,ISYM(MXAO),JCODE(MXAO),
     *                NLCS(MXAO),LEVPT(MXAO),LEVNR(MXAO),IOUT(MXAO),
     *                NREFS,IEXCT,NFOCI,INTACT
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FEX,UNZOR,PORT
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMEFP/ PB_MUL,PC_MUL,PB_POL,PC_POL,P_FF,P_NF,P_FN,P_I,
     *                P_IBIS,P_J,P_NUCC,P_NUCCBIS,ENPCM
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IKREP,IEF,IP_F,NFMOPCM
      COMMON /PCMPNT/ LXYZRE,LSSFE,LLIST,LALPCM,LRINPCM,LINAPCM,LINFPCM,
     *                LINIPCM,LMEPCM,LAXYZCT,LXYZCT2,LNVERT,LQSN,LQSE,
     *                LQFS,LQIND,LISPHE,LVAD,LQOR,LVECMUL,LAIPRJ,LFIPRJ
     *                ,LPEL
      COMMON /PCMPRT/ GCAVP,GCAVS,GDISP,GREP,EHFGAS
      COMMON /QUDMC / QUDTHR,QUD
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /SVPSCM/ IADDR_SVP,ILAST_SVP
      COMMON /SVPNRG/ ATMSVP,ELCSVP,SVPSVP,SURSVP,VOLPLM,VOLSVP,RXNSVP
      COMMON /SVPOPT/ ISVP,NSVP
      COMMON /SVPTGR/ NVLPL,MTHSVP,IQORP,NSTAGE,IOPPRD,IRENRM,ISHAPE,
     *                ITRCAV,NDRCAV,ICVICK,IFRQUP,IRDRF,IONSTR,
     *                INTCAV,NPTLEB,NPTTHE,NPTPHI,IROTGR,ITRNGR,
     *                ICRD,LINEQ,IPNRF,IPRRF,IPNDN,IPRDN,IEX
      COMMON /SVPPAR/ CAVC(3,MXCAV),CAVW(MXCAV),CAVQ(MXCAV),
     1                CAVN(3,MXCAV),SPHW(MXCAV),CAVP(MXCAV),
     2                CENCAV(3),RNTCAV(3),NCAV,IQP
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
      PARAMETER (CVGTOL=5.0D-03, DMPTOL=0.2D+00, TEN=1.0D+01,
     *           TEN7=1.0D-07, PT2=0.2D+00, TWOPT2=2.2D+00,
     *           ZERO=0.0D+00, TWO=2.0D+00, HALF=0.5D+00)
C
      SAVE INITAL
C
      DATA DETWRD,DRTWRD,GENWRD,ORMWRD
     *     /8HDET     ,8HDRT     ,8HGEN     ,8HORMAS   /
      DATA DET1,DET2/8HDET1    ,8HDET2    /
      DATA DRT1,DRT2/8HDRT1    ,8HDRT2    /
      DATA ALDET,GUGA,GENCI,ORMAS,GMCCI
     *     /8HALDET   ,8HGUGA    ,8HGENCI   ,8HORMAS   ,8HGMCCI   /
      DATA CHECK/8HCHECK   /
      DATA INITAL/.TRUE./
      DATA OPTMZE/8HOPTIMIZE/
      DATA ENGY/8HENERGY  /
      DATA   AMEX/8HMEX     /
      DATA DETNEO/8HNEODET  /
C
C     ----- MAIN DRIVER FOR MCSCF ORBITAL OPTIMIZATION -----
C     THIS CODE WAS WRITTEN BY MICHEL DUPUIS OF IBM-KINGSTON.
C     IT WAS ADAPTED FROM MICHEL'S HONDO7 PROGRAM TO GAMESS
C     AS THE FILES MCSCF.SRC AND MCTWO.SRC BY MIKE SCHMIDT
C     AT NDSU IN OCTOBER OF 1988.
C
      IF(QUD) THEN
         CALL MCQUD
         RETURN
      END IF
C
      NOTRF=.FALSE.
      SOME=NPRINT.NE.-5 .AND. MASWRK
      NPSAVE = NPRINT
      FULLNR = .NOT.QUD  .AND.  .NOT.SOSCF  .AND.  .NOT.FOCAS
     *                   .AND.  .NOT.JACOBI
C
      IF(IPCM.EQ.1) ENUCR=ENUCR-UNZ
      IF(SOME) THEN
         WRITE(IW,9000)
         IF(IPCM.NE.1) THEN
            WRITE(IW,9010) ENUCR
         ELSE
            WRITE(IW,9111) ENUCR,UNZ,ENUCR+UNZ
         END IF
      END IF
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
      IREST=2
      EMC =ZERO
      EMC0=ZERO
      SQCDF=ZERO
      ITER=0
      MICRO=0
      DE =ZERO
      DEAVG=ZERO
      EPSLON=TEN7
      LCNT=1
      ITSO=0
      IROTMX=0
      JROTMX=0
      IZRFLS = IZRF.NE.0
C
      IF(DAMP.LE.ZERO) DAMP=ZERO
      IF(FOCAS) DAMP=ZERO
      IF(SOSCF) DAMP=ZERO
      DMPING=DAMP.GT.ZERO
      CVGING=.FALSE.
      CVGED =.FALSE.
      CVGSVP    = .FALSE.
      CVGSVPDN  = .FALSE.
      CVGSVPDN1 = .FALSE.
      CVGSVPNR  = .FALSE.
      MICCVG=.FALSE.
      FINAL =.FALSE.
      MCNOCI=.FALSE.
      DOFOCK=.TRUE.
      DOEXCH=.FALSE.
C
      TOLLAG = ACURCY
      TOLENG = ENGTOL
      IF(RUNTYP.EQ.OPTMZE  .AND.  NEVALS.GT.0.AND.NFG.EQ.0) THEN
         GRMS = DDOT(3*NAT,EGRAD,1,EGRAD,1)
         GRMS = SQRT(GRMS/(3*NAT))
                               CFACT =  1.0D+00
         IF(GRMS.GT.0.005D+00) CFACT =  2.5D+00
         IF(GRMS.GT.0.020D+00) CFACT = 10.0D+00
         IF(GRMS.GT.0.100D+00) CFACT = 20.0D+00
         TOLLAG = CFACT*TOLLAG
         TOLENG = CFACT*TOLENG
      END IF
C
C        FROZEN ORBITAL RUNS DO NOT LEAD TO A SYMMETRIC LAGRANGIAN
C        SO THE CONVERGENCE TEST SHOULD ONLY BE ON THE ENERGY.
C
      IF(FCORE  .OR.  MOFRZ(1).NE.0  .OR.  NOROT(1,1).NE.0) THEN
         TOLNG2=TOLENG*TEN*TEN
         TOLLG2=TEN*TEN
      ELSE
         TOLNG2=TOLENG*TEN*TEN
         TOLLG2=TOLLAG*TWO
      END IF
C
      TIM0=ZERO
      TIM1=ZERO
      CALL TSECND(TIM0)
      TLEFTS = TIMLIM-TIM0
C
C     ----- DEFINE THE DETERMINANTS OR CONFIGURATIONS TO BE USED -----
C
      IF(NEORUN.EQ.1) THEN
         CALL EDETIN(NPFLG(1),DETWRD)
         CALL NDETIN(NPFLG(1),DETNEO)
         GO TO 50
      END IF
C
      IF(RUNTYP.NE.AMEX) THEN
         IF(CISTEP.EQ.ALDET) CALL DETINP(NPFLG(1),DETWRD)
         IF(CISTEP.EQ.GUGA)  CALL DRTGEN(NPFLG(1),DRTWRD)
         IF(CISTEP.EQ.GENCI) CALL GCIINP(NPFLG(1),GENWRD)
         IF(CISTEP.EQ.ORMAS) CALL FCINPUT(NPFLG(1),DETWRD,ORMWRD)
         IF(CISTEP.EQ.GMCCI) CALL OS_GMCCI_START(.TRUE.)
      ELSE
         IF(CISTEP.EQ.ALDET) THEN
            IF(MEXSTATE.EQ.1) CALL DETINP(NPFLG(1),DET1)
            IF(MEXSTATE.EQ.2) CALL DETINP(NPFLG(1),DET2)
         END IF
         IF(CISTEP.EQ.GUGA) THEN
            IF(MEXSTATE.EQ.1) CALL DRTGEN(NPFLG(1),DRT1)
            IF(MEXSTATE.EQ.2) CALL DRTGEN(NPFLG(1),DRT2)
         END IF
      END IF
   50 CONTINUE
C
C     AT PRESENT ONLY FOCAS AND SOSCF TAKE SOME ADVANTAGE OF THIS INFO.
      IF(NFG.NE.0.AND.NORBPROJ.NE.0) THEN
         IF(MASWRK) WRITE(IW,9024) NORBPROJ
      END IF
C
      IF (NORB.EQ.NUM) NORB = NQMT
      IF(CISTEP.EQ.ALDET .OR.
     *   CISTEP.EQ.GENCI .OR.
     *   CISTEP.EQ.ORMAS .OR.
     *   CISTEP.EQ.GMCCI) THEN
         MORB = NORB
         MORBS = NACT
         MCORBS = NCORSV
         IF(FULLNR) CALL SEQOPN(NFT16,'WORK16', 'UNKNOWN',
     *                          .FALSE.,'UNFORMATTED')
      END IF
      IF(CISTEP.EQ.GUGA) THEN
         MORB = NORB
         MORBS = NDOC+NALP+NAOS+NBOS+NVAL
         MCORBS = NMCC
      END IF
C
      IF (.NOT.(FOCAS.OR.SOSCF)  .AND.  MCORBS.EQ.0) THEN
         IF(SOME) WRITE(IW,*)
     *         'FORCING -DROPC- FALSE, THERE ARE NO CORE ORBITALS'
         DROPC=.FALSE.
      END IF
C         PARALLEL TRANSFORMATION WILL DIE IF THERE ARE NO VIRTUALS
C         SO ANY SUCH MBS + FVS JOB SHOULD BE RUN SERIALLY.
      IF(GOPARR  .AND.  FULLNR  .AND.  (MORB-MORBS-MCORBS).EQ.0) THEN
         IF(MASWRK) WRITE(IW,9026) MCORBS,MORBS,MORB
         CALL ABRT
      END IF
      IF (DROPC) MORBS = MORBS + MCORBS
      IF (DROPC) NORBS = NORBS + NCORBS
C
C         AT PRESENT, THE -DROPC- OPTION TO OMIT CORE ORBITALS IS
C         PROGRAMMED ONLY IN THE CODES THAT DON'T USE DISTRIBUTED
C         MEMORY TRANSFORMATIONS.  HENCE THE VARIABLE BELOW, WHICH
C         APPLIES ONLY TO THE DISTRIBUTED MEMORY TRANSFORM, JUST
C         KEEPS CORES IN THE TRANSFORM.
C
      DOCORE=.TRUE.
C
C        THE SETTING HERE PROBABLY APPLIES MAINLY TO THE ORBITAL
C        UPDATE STEPS, AS THE INTEGRAL TRANSFORMATION AND CI+DENSITY
C        STEPS TEND TO MANIPULATE DSKWRK AS NECESSARY.  IT'S TRICKY!
C
C        DURING ITERATIONS, GUGA HAS DUPLICATED FILES ON EACH NODE
C        BUT MOST PROCEDURES HAVE ONLY FILES ON THE MASTER NODE.
C
C        AFTER CONVERGENCE, DSKWRK SHOULD BE SET ON, SO THAT THE
C        FOCK MATRIX FOR CANONICALIZATION BUILDS CORRECTLY.
C
      DSKWRK = .FALSE.
C
C        THE FULL NEWTON-RAPHSON CONVERGER USES A DDI TRANSFORMATION
C
      IF (GOPARR  .AND.  FULLNR) THEN
         DDITRF = .TRUE.
         DOOOOO = .TRUE.
         DOVOOO = .TRUE.
         DOVVOO = .TRUE.
         DOVOVO = .TRUE.
         DOVVVO = .FALSE.
         DOVVVV = .FALSE.
      ELSE
         DDITRF = .FALSE.
         DOOOOO = .FALSE.
         DOVOOO = .FALSE.
         DOVVOO = .FALSE.
         DOVOVO = .FALSE.
         DOVVVO = .FALSE.
         DOVVVV = .FALSE.
      END IF
C
C     ----- OPEN -DIIS- WORK FILE FOR 1ST ORDER MCSCF -----
C
      IF (FOCAS .OR. SOSCF) THEN
         NCALL=0
         NDIIS=0
         L1=NUM
         L2 = (L1*L1+L1)/2
         L3 = L1*L1
         IDAF30=30
         NAV30 =0
         NDAR30=4*MDIIS+5
         LDAR30=MAX(L2,L3,MDIIS*MDIIS)
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL RAOPEN(IDAF30,IODA30,0,NDAR30,LDAR30,NPFLG(8))
         DSKWRK=SVDSKW
         DO 5 I=1, 2*MDIIS
            IDARE(I) = 0
            IDARF(I) = 0
    5    CONTINUE
      END IF
C
C     ----- OPEN DISK FILE FOR APPROXIMATE 2ND ORDER MCSCF -----
C
      IF(SOSCF) THEN
         IF(INITAL) THEN
            CALL SEQOPN(23,'JKFILE','UNKNOWN',.FALSE.,'UNFORMATTED')
            INITAL = .FALSE.
         ELSE
            CALL SEQREW(23)
         END IF
      END IF
C
C     ----- OPEN DISK FILE FOR JACOBI MCSCF-----
C
      IF(JACOBI) THEN
         IF(INITAL) THEN
            CALL SEQOPN(23,'JKFILE','UNKNOWN',.FALSE.,'UNFORMATTED')
         ELSE
            CALL SEQREW(23)
         END IF
      END IF
C
C     ----- EXCHANGE INTEGRALS FOR DIRECT TRANSFORMATION -----
C
      IF(DIRTRF) THEN
         NSH2 = (NSHELL*NSHELL+NSHELL)/2
         CALL BASCHK(LMAX)
                       NANGM =  4
         IF(LMAX.EQ.2) NANGM =  6
         IF(LMAX.EQ.3) NANGM = 10
         IF(LMAX.EQ.4) NANGM = 15
         IF(LMAX.EQ.5) NANGM = 21
         IF(LMAX.EQ.6) NANGM = 28
         MAXG = NANGM**4
         CALL VALFM(LOADFM)
         LXINTS = LOADFM + 1
         LGHOND = LXINTS + NSH2
         LDDIJ  = LGHOND + MAXG
         LAST   = LDDIJ  + 49*MXG2
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
         IF(SOME) WRITE(IW,9028)
         CALL EXCHNG(X(LXINTS),X(LGHOND),X(LDDIJ),NSH2,MAXG,INTTYP)
         CALL RETFM(NEED)
      END IF
C
C     -- SET UP COMBINED EFP/PCM
C
      IF(IPCM.EQ.1 .AND. IEFP.EQ.1) THEN
         IPCM=0
         IEFP=0
         IPCFP=.TRUE.
      ELSE
         IPCFP=.FALSE.
      END IF
C
C     -- MEMORY FOR SOLVATION METHODS
C
      CALL VALFM(LOADFM)
      LAST = LOADFM + 1
C
C     -- MEMORY FOR ANY OF THE SOLVATION METHODS
C
      IF (IZRFLS .OR. IEFP.EQ.1 .OR. IPCFP .OR.
     *    IPCM.EQ.1 .OR. NBUFMO.GT.0  .OR.  ISVP.EQ.1) THEN
         IF(ISVP.EQ.1) THEN
            LENW2=0
            LENW3=0
         ELSE
            LENW2=L2
            LENW3=L2
         END IF
         LD     = LAST
         LH1    = LD     + L2
         LC     = LH1    + L2
         LWRK1  = LC     + L3
         LWRK2  = LWRK1  + L2
         LWRK3  = LWRK2  + LENW2
         LAST   = LWRK3  + LENW3
      ELSE
         LD     = LAST
         LH1    = LAST
         LC     = LAST
         LWRK1  = LAST
         LWRK2  = LAST
         LWRK3  = LAST
      END IF
C
C      --- PCM AND PCM/EFP TOGETHER,
C        THE ORDER MUST BE PCM, EFP, PCM/EFP.
C
      CALL PCMMEM3(IPCFP,L2,LAST,NTSPAR,LSOL,LCHG,LDMATM,LQPOT,LQFLD,
     *             LVPOT,LSE,LDE,LQ0,LQ1,LQ2,LQ3,LD0,LQA,LDIMAT,LQREP,
     *             LTMP0,LTMP1,LTMP2,LIPVT,LPOTTMP,LRMUL,LCHG2,
     *             LASCCRD,LASCCHG,LASCDIP,LASCQAD,LDISV,LDIS1,LDIS2,
     *             LCQEF,LQEFF,LELD,LXDINT,LYDINT,LZDINT,LEFLD,MADD,
     *             LABFLD,LVECMP,LVEC_2,LVEC_S,LAPROJ,LFPROJ)
      LSCR   = LAST
      IF(IPCM.EQ.1.OR.IPCFP) THEN
         LCOL1  = LSCR   + L2
         LPOT   = LCOL1  + L1
         LBK    = LPOT   + L2
         LAST   = LBK    + L3
      ELSE
         LCOL1  = LAST
         LPOT   = LAST
         LBK    = LAST
      END IF
C
      IF(NFMOPCM.NE.0) THEN
         LDFMO  = LAST
         LH1FMO = LDFMO  + L2
         LAST   = LH1FMO + L2
      END IF
C
C         MEMORY FOR SVP SOLVENT MODEL
C
      IF(ISVP.EQ.1) THEN
         NOPER=2
         IF(IQP.EQ.0.AND.MTHSVP.EQ.1) NOPER=1
         LENSVP=10*NCAV+NOPER*NCAV**2
         LSVP=LAST
         LAST=LSVP+LENSVP
      ELSE
         LSVP=LAST
         LENSVP=0
      END IF
C          AND INITIALIZATIONS FOR -SVP- MODEL
      IUPRXN = -1
      IF(ISVP.EQ.1) THEN
         IADDR_SVP=LSVP
         ILAST_SVP=LSVP+LENSVP-1
         CALL SVPSET
      END IF
C
      IF (IZRFLS .OR. IEFP.EQ.1 .OR. IPCFP .OR. IPCM.EQ.1 .OR.
     *    NBUFMO.GT.0  .OR.  ISVP.EQ.1 .OR. NFMOPCM.NE.0) THEN
         NEED = LAST-LOADFM-1
         IF(MASWRK) WRITE(IW,9020) NEED
         CALL GETFM(NEED)
      END IF
C
C     -- READ INVERSE D MATRIX FOR PCM
C
      IF((IPCM.EQ.1.OR.IPCFP).AND.IPCMIT.EQ.0) THEN
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         NTS3 = NTS*NTS
         CALL SEQREW(NFT26)
         CALL DERCHK(NDER)
         IF(NDER.EQ.1  .AND.  IPCDER.EQ.1) THEN
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
         END IF
         IF(IEF.GT.0  .AND.  IEF.LT.3) THEN
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
            CALL SQREAD(NFT26,X(LSE),NTS3)
            CALL SQREAD(NFT26,X(LDE),NTS3)
         END IF
         CALL SQREAD(NFT26,X(LDMATM),NTS3)
         CALL SEQREW(NFT26)
         DSKWRK = SVDSKW
      END IF
C
C     ----- IF SCRF OR EFP OR PCM, SAVE PRISTINE ONE ELEC. HAM. -----
C     ----- IF SCRF READ IN DIPOLE INTEGRALS -----
C
      IF(IZRFLS .OR. IEFP.EQ.1. OR. IPCM.EQ.1 .OR.
     *   IPCFP .OR. NBUFMO.GT.0  .OR.  ISVP.EQ.1) THEN
         CALL DAREAD(IDAF,IODA,X(LH1),L2,11,0)
         IF(.NOT.(IPCFP.AND.RUNTYP.EQ.ENGY)) THEN
            CALL DAWRIT(IDAF,IODA,X(LH1),L2,87,0)
         END IF
      END IF
      IF(IZRFLS) THEN
         CALL DAREAD(IDAF,IODA,X(LXDINT),L2,95,0)
         CALL DAREAD(IDAF,IODA,X(LYDINT),L2,96,0)
         CALL DAREAD(IDAF,IODA,X(LZDINT),L2,97,0)
      END IF
C
C           **********************
C     ***** BEGIN MCSCF ITERATIONS *****
C           **********************
C
      IF(EXETYP.EQ.CHECK) MAXIT=1
  100 CONTINUE
      ITER=ITER+1
C
C     * * * * * SOLVATION MODEL CORRECTIONS TO 1E- TERM * * * * *
C
      IF (IZRFLS .OR. IEFP.EQ.1 .OR. IPCM.EQ.1 .OR.
     *    IPCFP .OR. NBUFMO.GT.0  .OR.  ISVP.EQ.1) THEN
         MCDENAO=.TRUE.
      ELSE
         MCDENAO=.FALSE.
      END IF
      IF(NBUFMO.GT.0) CALL DAREAD(IDAF,IODA,X(LH1),L2,87,0)
C
C                ----- ONSAGER REACTION FIELD -----
C
C     -- GET CURRENT DENSITY MATRIX AND PRISTINE ONE          --
C     -- ELEC. HAM., CALC. INDUCED DIPOLES, ADD PERTURBATION  --
C     --  LD   HOLDS TOTAL DENSITY,                           --
C     --  LH1  HOLDS CORE HAMILTONIAN PLUS SCRF PERTURBATION, --
C     -- AND IS WRITTEN TO DAF 11 INSIDE RCNFLD               --
C
      IF (IZRF.EQ.1.AND.ITER.GT.1) THEN
         CALL DAREAD(IDAF,IODA,X(LD),L2,16,0)
         CALL DAREAD(IDAF,IODA,X(LH1),L2,87,0)
         CALL RCNFLD(X(LH1),DUMMY,DUMMY,X(LD),DUMMY,DUMMY,
     *               X(LXDINT),X(LYDINT),X(LZDINT),L2,L1,ITER)
      END IF
C
C     ----- SURFACE AND VOLUME POLARIZATION FOR ELECTROSTATICS -----
C     ON THE FIRST ITERATION, WE DO NOT YET KNOW A DENSITY MATRIX
C     SO WE HAVE TO SKIP OVER THE -SVP- SOLVATION UNTIL 2ND ITER.
C
      IF(ISVP.EQ.1  .AND.  ITER.GT.1) THEN
         ITERSVP=ITER-1
         CALL DAREAD(IDAF,IODA,X(LD),L2,16,0)
         CALL DAREAD(IDAF,IODA,X(LH1),L2,87,0)
         CALL SVPCHG(X(LD),DUMMY,X(LH1),X(LWRK1),IADDR_SVP,ILAST_SVP,
     1               L2,L1,ESVP,ITERSVP,
     2               CVGSVPDN,CVGSVPDN1,CVGSVPNR,IUPRXN,CHGDIF,CVGSVP)
         CALL DAWRIT(IDAF,IODA,X(LH1),L2,11,0)
      END IF
C
C           ----- EFFECTIVE FRAGMENT POTENTIAL MODEL -----
C
C     -- GET CURRENT DENSITY MATRIX AND PRISTINE ONE          --
C     -- ELEC. HAM., CALC. INDUCED DIPOLES, ADD PERTURBATION  --
C     --  LD   HOLDS TOTAL DENSITY,                           --
C     --  LH1  HOLDS CORE HAMILTONIAN PLUS IND. DIPOLE        --
C     -- PERTURBATION, AND IS WRITTEN TO DAF 11 INSIDE EFPCMP --
C
      IF (IEFP.EQ.1.AND.ITER.GT.1) THEN
         CALL DAREAD(IDAF,IODA,X(LD),L2,16,0)
         CALL DAREAD(IDAF,IODA,X(LH1),L2,87,0)
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         IF(NBUFMO.GT.0) THEN
           CALL MCAODBF(X(LC),X(LWRK1),X(LWRK2),L1,NBUFMO)
           CALL VSUB(X(LWRK1),1,X(LD),1,X(LD),1,L2)
         END IF
         CALL EFPCMP(X(LH1),DUMMY,DUMMY,X(LD),DUMMY,DUMMY,
     *               X(LXDINT),X(LYDINT),X(LZDINT),X(LEFLD),
     *               X(MADD),X(LABFLD),L2,L1,ITER)
         DSKWRK = SVDSKW
      END IF
C
C     ----- POLARIZABLE CONTINUUM MODEL CALCULATION -----
C
C     -- GET CURRENT DENSITY MATRIX AND PRISTINE ONE          --
C     -- ELEC. HAM., ADD PERTURBATION                         --
C     --  LD   HOLDS TOTAL DENSITY,                           --
C     --  LH1  HOLDS CORE HAMILTONIAN PLUS PCM CHARGE         --
C     -- PERTURBATION, AND IS WRITTEN TO DAF 11 INSIDE PCMFLD --
C
      IF(IPCM.EQ.1.AND.ITER.GT.1) THEN
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
C
C        -- ADD PCM 1-E INTERACTION MATRIX
C
         IF(ITER.EQ.2.AND.IPCMIT.EQ.0) THEN
            IF(IP_F.EQ.1) THEN
             CALL VALFM(LOADFM)
             LVEC    = LOADFM   + 1
             LVECI   = LVEC     + L2
             LVECTEST= LVECI    + L2
             LAST    = LVECTEST + L2
             NEEDEF= LAST - LOADFM -1
             CALL GETFM(NEEDEF)
             CALL ENUADD(X(LDMATM),X(LVEC),X(LVECI),X(LVECTEST),
     *                   L2,X(LAXYZCT),
     *                   X(LAXYZCT+MXTS),X(LAXYZCT+MXTS*2),
     *                   X(LAXYZCT+MXTS*3),X(LQSN),X(LQFS),X(LVECMUL))
             CALL DAWRIT(IDAF,IODA,X(LVECI),L2,263,0)
             CALL RETFM(NEEDEF)
            END IF
C
            ENUCR=ENUCR+UNZ
            CALL JMAT(X(LSOL),X(LSCR),X(LPOT),X(LDIS1),X(LDIS2),L2,
     *                X(LAXYZCT),X(LAXYZCT+MXTS),X(LAXYZCT+MXTS*2),
     *                X(LQSN),X(LQOR))
            CALL DAWRIT(IDAF,IODA,X(LSCR),L2,264,0)
         END IF
C
         CALL DAREAD(IDAF,IODA,X(LD) ,L2,16,0)
         CALL DAREAD(IDAF,IODA,X(LH1),L2,87,0)
         CALL DCOPY(L2,X(LD),1,X(LWRK1),1)
         IF(NBUFMO.GT.0) THEN
           CALL MCAODBF(X(LC),X(LWRK1),X(LWRK2),L1,NBUFMO)
           CALL VSUB(X(LWRK1),1,X(LD),1,X(LWRK1),1,L2)
         END IF
         CALL PCMFLD(X(LH1),DUMMY,X(LD),X(LSOL),X(LCHG),X(LSCR),X(LDISV)
     *              ,X(LDIS1),X(LDIS2),X(LBK),X(LCOL1),X(LDMATM),
     *               X(LCQEF),X(LQEFF),X(LELD),X(LQPOT),X(LQFLD),
     *               X(LVPOT),X(LSE),X(LDE),
     *               X(LQ0),X(LQ1),X(LQ2),X(LQ3),
     *               X(LD0),X(LQA),X(LDIMAT),X(LQREP),
     *               X(LTMP0),X(LTMP1),X(LTMP2),X(LIPVT),
     *               X(LPOTTMP),X(LRMUL),TCH,QET,QETN,QESC,1,NFT27,
     *               .FALSE.,L1,L2,NTSPAR,X(LAXYZCT),X(LAXYZCT+MXTS),
     *               X(LAXYZCT+MXTS*2),X(LXYZCT2),X(LXYZCT2+MXTS),
     *               X(LXYZCT2+MXTS*2),X(LAXYZCT+MXTS*3),
     *               X(LQSN),X(LQSE),X(LQFS),X(LVECMUL),
     *               X(LCHG2),X(LQIND+MXTS),X(LPEL),X(LASCCRD),
     *               X(LASCCHG),X(LASCDIP),X(LASCQAD),X(LXYZRE),
     *               X(LXYZRE+MXSP),X(LXYZRE+MXSP*2),X(LXYZRE+MXSP*3),
     *               X(LISPHE),X(LLIST))
         DSKWRK = SVDSKW
      END IF
C
C               ------ COMBINED PCM/EFP METHOD ------
C
C     -- GET CURRENT DENSITY MATRIX AND PRISTINE               --
C     -- ONE ELEC. HAM., ADD PERTURBATION                      --
C     --  LD   HOLDS TOTAL DENSITY,                            --
C     --  LH1  HOLDS CORE HAMILTONIAN PLUS CHARGE AND DIPOLE   --
C     -- PERTURBATION, AND IS WRITTEN TO 11 IN PCMFLD & CHGIND --
C
      IF(IPCFP.AND.ITER.GT.1) THEN
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
C
C        -- ADD PCM 1-E INTERACTION MATRIX
C
         IF(ITER.EQ.2.AND.IPCMIT.EQ.0) THEN
            IF(IP_F.EQ.1) THEN
             CALL VALFM(LOADFM)
             LVEC    = LOADFM  + 1
             LVECI   = LVEC    + L2
             LVECTEST= LVECI   + L2
             LAST    = LVECTEST+ L2
             NEEDEF = LAST - LOADFM -1
             CALL GETFM(NEEDEF)
             CALL ENUADD(X(LDMATM),X(LVEC),X(LVECI),X(LVECTEST),
     *                   L2,X(LAXYZCT),
     *                   X(LAXYZCT+MXTS),X(LAXYZCT+MXTS*2),
     *                   X(LAXYZCT+MXTS*3),X(LQSN),X(LQFS),X(LVECMUL))
             CALL DAWRIT(IDAF,IODA,X(LVECI),L2,263,0)
             CALL RETFM(NEEDEF)
            END IF
C
            ENUCR=ENUCR+UNZ
            CALL JMAT(X(LSOL),X(LSCR),X(LPOT),X(LDIS1),X(LDIS2),L2,
     *                X(LAXYZCT),X(LAXYZCT+MXTS),X(LAXYZCT+MXTS*2),
     *                X(LQSN),X(LQOR))
            CALL DAWRIT(IDAF,IODA,X(LSCR),L2,264,0)
         END IF
C
         CALL DAREAD(IDAF,IODA,X(LD) ,L2,16,0)
         CALL DAREAD(IDAF,IODA,X(LH1),L2,87,0)
C
         CALL DCOPY(L2,X(LD),1,X(LWRK1),1)
         IF(NBUFMO.GT.0) THEN
           CALL MCAODBF(X(LC),X(LWRK1),X(LWRK2),L1,NBUFMO)
           CALL VSUB(X(LWRK1),1,X(LD),1,X(LWRK1),1,L2)
         END IF
C
         CALL PCMFLD(X(LH1),DUMMY,X(LD),X(LSOL),X(LCHG),X(LSCR),X(LDISV)
     *              ,X(LDIS1),X(LDIS2),X(LBK),X(LCOL1),X(LDMATM),
     *               X(LCQEF),X(LQEFF),X(LELD),X(LQPOT),X(LQFLD),
     *               X(LVPOT),X(LSE),X(LDE),
     *               X(LQ0),X(LQ1),X(LQ2),X(LQ3),
     *               X(LD0),X(LQA),X(LDIMAT),X(LQREP),
     *               X(LTMP0),X(LTMP1),X(LTMP2),X(LIPVT),
     *               X(LPOTTMP),X(LRMUL),TCH,QET,QETN,QESC,1,NFT27,
     *               IPCFP,L1,L2,NTSPAR,X(LAXYZCT),X(LAXYZCT+MXTS),
     *               X(LAXYZCT+MXTS*2),X(LXYZCT2),X(LXYZCT2+MXTS),
     *               X(LXYZCT2+MXTS*2),X(LAXYZCT+MXTS*3),
     *               X(LQSN),X(LQSE),X(LQFS),X(LVECMUL),
     *               X(LCHG2),X(LQIND+MXTS),X(LPEL),X(LASCCRD),
     *               X(LASCCHG),X(LASCDIP),X(LASCQAD),X(LXYZRE),
     *               X(LXYZRE+MXSP),X(LXYZRE+MXSP*2),X(LXYZRE+MXSP*3),
     *               X(LISPHE),X(LLIST))
C
         CALL CHGIND(X(LH1),DUMMY,DUMMY,X(LWRK1),DUMMY,DUMMY,
     *               X(LDMATM),X(LXDINT),X(LYDINT),X(LZDINT),
     *               X(LEFLD),X(MADD),X(LABFLD),X(LVECMP),X(LVEC_2),
     *               X(LVEC_S),L2,L1,ITER,
     *               X(LAXYZCT),X(LAXYZCT+MXTS),X(LAXYZCT+MXTS*2),
     *               X(LAXYZCT+MXTS*3),X(LQSN),X(LQSE),X(LQFS),
     *               X(LQIND),X(LAIPRJ),X(LFIPRJ))
         DSKWRK = SVDSKW
      END IF
C
C
C   SUBTRACT INTERACTION OF THE EFP WITH THE BUFFER ZONE DENSITY.
C
C   HERE WE DEFINE EEFP AS THE INTERACTION OF THE EFP(MULTIPOLE
C   POINTS AND POLARIZABLE POINTS) WITH THE BUFFER ZONE DENSITY
C   I.E. (1)BUFFER ELECTRON -- MULTIPOLE POINTS
C        (2)BUFFER ELECTRON -- POLARIZABLE POINTS(INDUCED DIPOLE)
C   ANY PCM ASC ASSOCIATED INTERACTION IS NOT CONSIDERED AS A
C   PART OF EEFP.
C
      IF(NBUFMO.GT.0.AND.ITER.GT.1) THEN
        CALL DAREAD(IDAF,IODA,X(LWRK1),L2,319,0)
        CALL VSUB(X(LWRK1),1,X(LH1),1,X(LWRK1),1,L2)
        IF(IPCFP.OR.IPCM.EQ.1)
     *     CALL VSUB(X(LSOL),1,X(LWRK1),1,X(LWRK1),1,L2)
C
        IF(IPCMIT.EQ.0) THEN
           IF(IPCM.EQ.1.OR.IPCFP) THEN
              CALL DAREAD(IDAF,IODA,X(LSCR),L2,263,0)
              CALL VSUB(X(LSCR),1,X(LWRK1),1,X(LWRK1),1,L2)
              CALL DAREAD(IDAF,IODA,X(LSCR),L2,264,0)
              CALL VSUB(X(LSCR),1,X(LWRK1),1,X(LWRK1),1,L2)
           END IF
           IF(IPCFP) THEN
              CALL VSUB(X(LVEC_2),1,X(LWRK1),1,X(LWRK1),1,L2)
           END IF
        END IF
C
        CALL MCAODBF(X(LC),X(LWRK3),X(LWRK2),L1,NBUFMO)
        EEFP= TRACEP(X(LWRK3),X(LWRK1),L1)
      ELSE
        EEFP = ZERO
      END IF
      IF(NFMOPCM.NE.0) THEN
        CALL DAREAD(IDAF,IODA,X(LDFMO) ,L2,16,0)
        CALL DAREAD(IDAF,IODA,X(LH1FMO),L2,311,0)
        PX=-TRACEP(X(LDFMO),X(LH1FMO),L1)
      ENDIF
C
C     ----- GET CI WAVEFUNCTION AND ENERGY (E IS SET BY GUG2DM) -----
C
C        FIRST, PERFORM APPROPRIATE INTEGRAL TRANSFORMATION
C
                                SVDSKW=DSKWRK
                                DSKWRK=.TRUE.
      IF(FULLNR  .AND.  GOPARR) DSKWRK=.FALSE.
C        PURIFICATION IS DONE IN TRFMCX IF AND ONLY IF IPURTF=2,
C        SO FOR "FIRST ITERATION ONLY" PURIFICATION, RESET THIS.
      IPURSAVE = IPURTF
      IF(IPURTF.EQ.1  .AND.  ITER.EQ.1) IPURTF = 2
C
      IF (FOCAS .OR. SOSCF) THEN
         CALL CASTRF(NPFLG(2),.FALSE.)
      END IF
      IF(JACOBI) THEN
         CALL TRFMCX(NPFLG(2),MCORBS,MORBS,MORBS,DOFOCK,DOEXCH,
     *               DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *               DOVVVO,DOVVVV,DOCORE)
      END IF
      IF(FULLNR) THEN
         IF (DROPC) THEN
            CALL TRFMCX(NPFLG(2),MCORBS,MORBS,MORB,DOFOCK,DOEXCH,
     *                  DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *                  DOVVVO,DOVVVV,DOCORE)
         ELSE
            CALL TRFMCX(NPFLG(2),0,0,0,.FALSE.,.FALSE.,
     *                  DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *                  DOVVVO,DOVVVV,DOCORE)
         END IF
      END IF
      DSKWRK=SVDSKW
      IPURTF=IPURSAVE
C
C        EVALUATE CI WAVEFUNCTION AND GENERATE TWO PARTICLE DENSITY
C
      IF(NEORUN.EQ.1) THEN
         CALL NEOFCI(NPFLG,ITER-1)
         GO TO 60
      END IF
C
      IF(CISTEP.EQ.GUGA) THEN
         SVDSKW=DSKWRK
         DSKWRK=.TRUE.
         CALL GUGSRT(NPFLG(3),DDITRF)
         CALL GUGAEM(NPFLG(3))
         CALL GUGADG(NPFLG(4))
         CALL GUG2DM(NPFLG(6))
         DSKWRK=SVDSKW
      END IF
      IF(CISTEP.EQ.ALDET) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL DETFCI(NPFLG(3),.TRUE.,DDITRF)
         CALL DETDM2(NPFLG(6))
         DSKWRK=SVDSKW
      END IF
      IF(CISTEP.EQ.GENCI) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL DETGCI(NPFLG(3),.TRUE.,DDITRF)
         CALL GCIDM2(NPFLG(6))
         DSKWRK=SVDSKW
      END IF
      IF(CISTEP.EQ.ORMAS) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL DEFCCI(NPFLG(3),.TRUE.,DDITRF,0,1,0,NPFLG(6))
         DSKWRK=SVDSKW
      END IF
      IF(CISTEP.EQ.GMCCI) THEN
         SVDSKW=DSKWRK
         DSKWRK=.FALSE.
         CALL OS_GMCCI_SOLVE(DDITRF)
         DSKWRK=SVDSKW
      END IF
   60 CONTINUE
C
C        OPTIONAL LINEAR SEARCH IN ORBITAL SPACE (FULLNR METHOD ONLY)
C
      IF(LINSER) CALL LINMCM(EMC,ITER,LCNT,DROPC,MORBS,MCORBS,
     *                       DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *                       DOVVVO,DOVVVV,DOCORE)
C
C        FINAL SOLVENT ENERGY CORRECTIONS
C     NOTE: ONSAGER MODEL PUT FREE ENERGY CORRECTIONS INTO -EN- ALREADY
C
      IF(ITER.GT.1) E = E - EEFP
C
      IF((IPCFP.OR.IPCM.EQ.1).AND.ITER.GT.1) THEN
        IF(IPCMIT.EQ.0) THEN
           E=E-0.50D+00*(PB+PX-PC)
        ELSE
           E=E+PX-0.50D+00*ENPCM
        END IF
      END IF
      IF(NFMOPCM.NE.0) E=E+0.5D+00*(PC+PX)
C
      IF(IPCM.EQ.1) THEN
        E=E+(GCAVP+GDISP+GREP)/627.509541D+00
      END IF
C
C     SVP MODEL'S FREE ENERGY CORRECTION IS
C     +HALF*ATMSVP SINCE IT HAS NOT YET BEEN ACCOUNTED FOR, BUT
C     -HALF*ELCSVP SINCE EHF ALREADY IMPLICITLY CONTAINS ELCSVP
C
      IF(ISVP.EQ.1) THEN
         E=E+HALF*(ATMSVP-ELCSVP)
      END IF
C
      EMC0=EMC
      EMC =E
      ETOT=E
      DE0=DE
      DE =EMC-EMC0
      IF(ITER.EQ.1) DEAVG=ZERO
      IF(ITER.EQ.2) DEAVG= ABS(DE)
      IF(ITER.GE.3) DEAVG=( ABS(DE)+ ABS(DE0)+PT2*DEAVG)/TWOPT2
C
      EPSLON=EPSLON/TEN
      IF(ABS(DE).LE.ACURCY) EPSLON=TEN7
      IF(EPSLON.LT.TEN7) EPSLON=TEN7
      TOLMIC= ABS(EMC)* SQRT(EPSLON)
C
      CVGING=ABS(DE).LT.CVGTOL .AND. DAMP.LT.DMPTOL
      IF(DMPING.AND.     CVGING                 ) DAMP=ZERO
      IF(DMPING.AND..NOT.CVGING.AND.DAMP.EQ.ZERO) DAMP=DMPTOL
C
      IF(.NOT.SOSCF .AND. .NOT.FOCAS) THEN
         IF(ITER.GT.2.AND..NOT.CVGING) CALL NTNDMP(DE,DE0,DEAVG,DAMP)
         IF(.NOT.DMPING.AND.DAMP.GT.ZERO) DMPING=.TRUE.
      END IF
C
C     ***** START MICRO-ITERATION *****
C
      ENERG0=ZERO
      MICRO=0
  500 CONTINUE
      MICRO=MICRO+1
C
C     ----- MCSCF ORBITAL IMPROVEMENT -----
C
      DEAMAX=ZERO
      IF(FOCAS .OR. SOSCF) THEN
         SVDSKW=DSKWRK
         DSKWRK=.TRUE.
         CALL CASMAC(NPFLG(8),.FALSE.,ITSO,IROTMX,JROTMX,DEAMAX)
         DSKWRK=SVDSKW
      END IF
      IF(JACOBI) THEN
         CALL JIKRJA(INITAL,NTROT,ECON)
         INITAL=.FALSE.
      END IF
      IF(FULLNR) THEN
         CALL NEWTON(CISTEP,NPFLG(8),METHOD,ACURCY,DAMP,EKT,
     *               CANONC,NPUNCH,FCORE,FORS,NOCI,
     *               NWORD,MORB,NOROT,MXNORO,IROTMX,JROTMX,MOFRZ,
     *               FINAL,MCNOCI,DOEXCH,DROPC,FULLNR)
      END IF
C
C     --- IF NEEDED TRANSFORM MCSCF TO AO BASIS AND WRITE TO DAF 16 ---
C
      IF (MCDENAO) CALL MCAODEN(X(LC),X(LD),X(LH1),L1,MORBS)
C
      CALL TSECND(TIM1)
      TIM0=TIM1
      IF(MASWRK  .AND. .NOT.JACOBI) THEN
         IF(ITER.EQ.1) WRITE(IW,9030)
         IF(ITSO.EQ.1) WRITE(IW,9035)
         IHART = INT(ABS(DE))
         IF(DE.LE.ZERO) THEN
            DELE=DE+IHART
         ELSE
            DELE=DE-IHART
         END IF
         WRITE(IW,9040) ITER,EMC,DELE,DEMAX,IROTMX,JROTMX,SQCDF,
     *                  MICRO,DAMP
      END IF
      CALL FLSHBF(IW)
C
C        NOW SUPPRESS FULL PRINTOUT ON FOLLOWING ITERATIONS.
C
      IF(NPFLG(10).EQ.0) THEN
         NPRINT = -5
         IF(NPFLG(1).LE.0) NPFLG(1)=-5
         IF(NPFLG(2).LE.0) NPFLG(2)=-5
         IF(NPFLG(3).LE.0) NPFLG(3)=-5
         IF(NPFLG(4).LE.0) NPFLG(4)=-5
         IF(NPFLG(6).LE.0) NPFLG(6)=-5
         IF(NPFLG(8).LE.0) NPFLG(8)=-5
      END IF
C
C   TEST FOR CONVERGENCE FOR JACOBI METHOD
C
      IF (JACOBI.AND.ECON) GO TO 1000
C
      IF(JACOBI) GO TO 1314
C
      MICCVG=DEMAX.LT.TOLMIC.OR.ENERGY.GT.ENERG0
      ENERG0=ENERGY
      IF(MICCVG.OR.MICRO.GE.MICIT) GO TO 510
      IF(FOCAS .OR. SOSCF) GO TO 510
                                SVDSKW=DSKWRK
                                DSKWRK=.TRUE.
      IF(FULLNR  .AND.  GOPARR) DSKWRK=.FALSE.
      IF (DROPC) THEN
         CALL TRFMCX(NPFLG(2),MCORBS,MORBS,MORB,DOFOCK,DOEXCH,
     *               DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *               DOVVVO,DOVVVV,DOCORE)
      ELSE
         CALL TRFMCX(NPFLG(2),0,0,0,.FALSE.,.FALSE.,
     *               DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *               DOVVVO,DOVVVV,DOCORE)
      END IF
      DSKWRK=SVDSKW
      GO TO 500
C
C     ***** END OF MICRO-ITERATIONS *****
C
  510 CONTINUE
C
C     ----- TEST FOR CONVERGENCE -----
C
      IF (ISTAT.EQ.1) GO TO 1314
C
      CVGED = CVGING .AND. ABS(DE).LT.TOLNG2 .AND. DEMAX.LT.TOLLAG
C
      IF(ISVP.EQ.1) THEN
         CVGED  = CVGED   .AND.  CVGSVP
         CVGING = CVGING  .AND.  CVGSVP
         CVGSVPNR = DEMAX  .LT.  10.0D+00 * TOLLAG
         CVGSVPDN = CVGING
         CVGSVPDN1 = CVGING
      END IF
C
      IF(CVGED .AND. MASWRK) WRITE(IW,9050)
      IF(CVGED) GO TO 1000
      CVGED = CVGING .AND. ABS(DE).LT.TOLENG .AND. DEMAX.LT.TOLLG2
      IF(CVGED .AND. MASWRK) WRITE(IW,9060)
      IF(CVGED) GO TO 1000
      CVGED = CVGING .AND.  ABS(DE).LT.TOLNG2
     *               .AND.  DEMAX.LT.5*TOLLAG
     *               .AND.  SQCDF.LT.1.0D-11
      IF(CVGED .AND. MASWRK) WRITE(IW,9065)
      IF(CVGED) GO TO 1000
C
C     ----- TEST FOR TIME LIMIT -----
C
 1314 CONTINUE
      CALL TSECND(TIM1)
      TLEFT = TIMLIM-TIM1
      IF(TLEFT .LT.  1.5D+00*(TLEFTS-TLEFT)/ITER) THEN
         IF (MASWRK) THEN
            WRITE(IW,9070)
            WRITE(IW,9090)
         END IF
         EMC=ZERO
         E  =ZERO
         GO TO 1000
      END IF
C
C     ---- TEST FOR ITERATION LIMIT -----
C
      IF(ITER.GE.MAXIT) THEN
         IF (MASWRK) THEN
            WRITE(IW,9080)
            WRITE(IW,9090)
         END IF
         EMC=ZERO
         E  =ZERO
         GO TO 1000
      END IF
      GO TO 100
C
C           **********************
C     ***** END OF MCSCF ITERATION *****
C           **********************
C
 1000 CONTINUE
C
      IF(NFG.NE.0.AND.MCFMO.NE.0.AND.DEAMAX.GT.DEMAX.AND.MASWRK) THEN
        IF(MCFMO.EQ.1) WRITE(IW,9075) DEAMAX
        IF(MCFMO.EQ.2) WRITE(IW,9077) DEAMAX
      ENDIF
C
      DAMP = 0.0D+00
C
C     IF THIS WAS A STATE AVERAGED MCSCF, THE USER MAY WANT THE ENERGY
C     SET TO SOME SPECIFIC STATE'S VALUE RATHER THAN THE AVERAGED ENERGY
C
      IF(IROOT.GT.0) THEN
         IF(CISTEP.EQ.ALDET  .OR.
     *      CISTEP.EQ.GENCI  .OR.
     *      CISTEP.EQ.ORMAS  .OR.
     *      CISTEP.EQ.GMCCI) THEN
            NXTR=0
            DO 310 ISTAT=1,KDET
               IF(IPURES.EQ.1  .AND.
     *            ABS(SPINS(ISTAT)-SDET).GT.0.03D+00) GO TO 310
               NXTR=NXTR+1
               IF(NXTR.EQ.IROOT) ESPEC=ESTATE(ISTAT)
  310       CONTINUE
         END IF
         IF(CISTEP.EQ.GUGA) THEN
            ESPEC = ESTATE(IROOT)
         END IF
         E = ESPEC
      END IF
C
C     ----- PRINT WAVEFUNCTION -----
C
      NPRINT = NPSAVE
      IF(EMC.NE.ZERO) EMC = E
      IF (MASWRK) WRITE(IW,9100) EMC,ITER
C
C         SOLVATION OUTPUT FOR PCM MODEL
C
      IF(IPCFP) THEN
         IPCM=1
         IEFP=1
         IPCFP=.FALSE.
      END IF
C
      QTT=0.0D+00
      DO I=0,NTS-1
C        QT=QSE(I)+QSN(I)+Q_FS(I)+Q_IND(I)
         QT=X(LQSE+I)+X(LQSN+I)+X(LQFS+I)+X(LQIND+I)
         QTT=QTT+QT
      ENDDO
      IF(IPCM.EQ.1.AND.IPCMIT.EQ.1) THEN
        IF(MASWRK) WRITE(IW,*) ' '
        IF(MASWRK) WRITE(IW,'(A10,F18.10)')'PCM ASC=',QTT
      END IF
      IF(IPCM.EQ.1) THEN
        IF(MASWRK) WRITE(IW,'(A10,F18.10)')
     *      'GCAVP  =',GCAVP/627.509541D+00
        IF(MASWRK) WRITE(IW,'(A10,F18.10)')
     *      'GDISP  =',GDISP/627.509541D+00
        IF(MASWRK) WRITE(IW,'(A10,F18.10)')
     *      'GREP   =',GREP/627.509541D+00
      END IF
C
      IF(MASWRK  .AND.  IPCM.EQ.1  .AND.  IPCMIT.EQ.0) THEN
         WRITE(IW,*)
         WRITE(IW,*)'-- INDUCED PCM ELECTRONIC CHARGES --'
         IF(ICOMP.EQ.0) WRITE(IW,1010) QET, TCH
         IF(ICOMP.EQ.2) WRITE(IW,1100) QET, TCH, QETN
         IF(ICOMP.EQ.3) WRITE(IW,1110) QET, TCH, QESC, QETN
         WRITE(IW,*)'--------------------------------------'
      END IF
C
C         SOLVATION OUTPUT FOR SVP MODEL
C
      IF(ISVP.EQ.1  .AND.  MASWRK) THEN
         CALL SVPOUT(ETOT,ECHIP,NPRINT,'MCSCF ')
      END IF
C
C         SOLVATION OUTPUT FOR ONSAGER OR EFP MODELS
C
      IF ((IZRF.EQ.1 .OR. IEFP.EQ.1) .AND.  MASWRK)
     *     WRITE(IW,9135) TMUX,TMUXD,TMUY,TMUYD,TMUZ,TMUZD
C
C     --- IF NEO, ONE LAST NEO-CI, OUPUT COEFFICIENTS, NO'S, RETURN ---
C
      IF(NEORUN.EQ.1) THEN
         CALL CASTRF(NPFLG(2),.FALSE.)
         CALL NEOFCI(NPFLG,-1)
C
         IF(CVGED) IREST=0
         IF (MASWRK) WRITE(IW,9110)
         CALL TIMIT(1)
C
         IF(FOCAS .OR. SOSCF) THEN
            SVDSKW = DSKWRK
            DSKWRK = .FALSE.
            CALL RACLOS(IDAF30,'DELETE')
            DSKWRK = SVDSKW
         END IF
         RETURN
      END IF
C
C         OUTPUT OF CONVERGED ORBITALS, INCLUDING CANONICALIZATION
C
      FINAL  = .TRUE.
      SVDSKW = DSKWRK
      DSKWRK = .TRUE.
      CALL NEWTON(CISTEP,NPFLG(8),METHOD,ACURCY,DAMP,EKT,
     *            CANONC,NPUNCH,FCORE,FORS,NOCI,
     *            NWORD,MORB,NOROT,MXNORO,IROTMX,JROTMX,MOFRZ,
     *            FINAL,MCNOCI,DOEXCH,DROPC,FULLNR)
      DSKWRK = SVDSKW
C
C     ----- END OF MCSCF -----
C
      IF(CVGED) IREST=0
      IF (MASWRK) WRITE(IW,9110)
      CALL TIMIT(1)
C
C     ----- CLOSE -DIIS- WORK FILE -----
C
      IF(FOCAS .OR. SOSCF) THEN
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         CALL RACLOS(IDAF30,'DELETE')
         DSKWRK = SVDSKW
      END IF
C
C     ----- DO NATURAL ORBITAL -CI- IF REQUESTED -----
C
      IF(NOCI  .AND.  CVGED) THEN
         IF (MASWRK) WRITE(IW,9120)
         SVDSKW=DSKWRK
         DSKWRK=.TRUE.
         IF(FULLNR  .AND.  GOPARR) DSKWRK=.FALSE.
         IF (FOCAS .OR. SOSCF) THEN
            CALL CASTRF(0,.FALSE.)
         ELSE
            IF (DROPC) THEN
               CALL TRFMCX(0,MCORBS,MORBS,MORB,DOFOCK,DOEXCH,
     *               DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *               DOVVVO,DOVVVV,DOCORE)
            ELSE
               CALL TRFMCX(0,0,0,0,.FALSE.,.FALSE.,
     *               DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *               DOVVVO,DOVVVV,DOCORE)
            END IF
         END IF
         DSKWRK=SVDSKW
C
         IF(CISTEP.EQ.GUGA) THEN
            SVDSKW=DSKWRK
            DSKWRK=.TRUE.
            CALL GUGSRT(0,DDITRF)
            CALL GUGAEM(0)
            CALL GUGADG(0)
            CALL GUG2DM(0)
            DSKWRK=SVDSKW
         END IF
         IF(CISTEP.EQ.ALDET) THEN
            SVDSKW=DSKWRK
            DSKWRK=.FALSE.
            CALL DETFCI(NPFLG(3),.TRUE.,DDITRF)
            CALL DETDM2(NPFLG(6))
            DSKWRK=SVDSKW
         END IF
         IF(CISTEP.EQ.GENCI) THEN
            SVDSKW=DSKWRK
            DSKWRK=.FALSE.
            CALL DETGCI(NPFLG(3),.TRUE.,DDITRF)
            CALL GCIDM2(NPFLG(6))
            DSKWRK=SVDSKW
         END IF
         IF(CISTEP.EQ.ORMAS) THEN
            SVDSKW=DSKWRK
            DSKWRK=.FALSE.
            CALL DEFCCI(NPFLG(3),.TRUE.,DDITRF,0,1,.TRUE.,-5)
            DSKWRK=SVDSKW
         END IF
         IF(CISTEP.EQ.GMCCI) THEN
            SVDSKW=DSKWRK
            DSKWRK=.FALSE.
            CALL OS_GMCCI_SOLVE(DDITRF)
            DSKWRK=SVDSKW
         END IF
         MCNOCI=.TRUE.
         CALL NEWTON(CISTEP,NPFLG(8),METHOD,ACURCY,DAMP,EKT,
     *               CANONC,NPUNCH,FCORE,FORS,NOCI,
     *               NWORD,MORB,NOROT,MXNORO,IROTMX,JROTMX,MOFRZ,
     *               FINAL,MCNOCI,DOEXCH,DROPC,FULLNR)
      END IF
C
C         ANALYTIC HESSIAN NEEDS TO REOPEN FILE 9 ON ALL NODES,
C         SO DELETE THE ONE ELECTRON MO INTEGRAL ON MASTER NODE.
C
      IF(GOPARR  .AND.  FULLNR) THEN
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         CALL SEQCLO(IJKT,'DELETE')
         DSKWRK = SVDSKW
      END IF
C
      IF (IZRFLS .OR. IEFP.EQ.1 .OR. IPCM.EQ.1 .OR. NBUFMO.GT.0 .OR.
     *    ISVP.EQ.1 .OR. NFMOPCM.NE.0) THEN
         CALL RETFM(NEED)
      END IF
      RETURN
C
 9000 FORMAT(/10X,17(1H-)/10X,'MCSCF CALCULATION'/10X,17(1H-))
 9010 FORMAT(/' ----- NUCLEAR ENERGY ----- = ',F20.10)
 9111 FORMAT(1X,'PCM SOLVATION MODEL:'/
     *       1X,'               NUCLEAR REPULSION ENERGY=',F20.10/
     *       1X,'NUCLEI-NUCLEAR POLARIZATION INTERACTION=',F20.10/
     *       1X,'                   TOTAL NUCLEAR ENERGY=',F20.10)
 9020 FORMAT(1X,'SOLVATION MODEL FOR MCSCF REQUIRES',I10,
     *    ' ADDITIONAL WORDS OF MEMORY')
 9024 FORMAT(/1X,'THE NUMBER OF PROJECTED OUT ORBITALS IN FMO IS',I4/)
 9026 FORMAT(1X,'ERROR, PARALLEL FULLNR RUN HAS NO VIRTUAL ORBITALS.'/
     *       1X,'PLEASE RUN IN SERIAL MODE TO DO THIS CORE/ACTIVE',
     *          ' OPTIMIZATION.'/
     *       1X,'NO.CORE=',I5,'  NO.ACT=',I5,' NO.TOTAL MO-S=',I5)
 9028 FORMAT(/1X,'EVALUATING EXCHANGE INTEGRALS FOR DIRECT',
     *           ' INTEGRAL TRANSFORMATION')
 9030 FORMAT(/1X,'ITER',5X,'TOTAL ENERGY',6X,'DEL(E)',
     *        2X,'LAGRANGIAN ASYMMETRY',3X,'SQCDF',2X,'MICIT',
     *        3X,'DAMP')
 9035 FORMAT(7X,10(1H-),'START APPROXIMATE SECOND ORDER MCSCF',10(1H-))
 9040 FORMAT(1X,I3,F19.9,F13.9,F10.6,2I5,1P,E10.3,0P,I3,F9.4)
 9050 FORMAT(/10X,20(1H-)/10X,'LAGRANGIAN CONVERGED'/10X,20(1H-))
 9060 FORMAT(/10X,16(1H-)/10X,'ENERGY CONVERGED'/10X,16(1H-))
 9065 FORMAT(/10X,15(1H-)/10X,'SQCDF CONVERGED (BUT THE LAGRANGIAN',
     *           ' MATRIX IS NOT ENTIRELY SYMMETRIC)'/10X,15(1H-))
 9070 FORMAT(1X,'RUNNING OUT OF CPU TIME...')
 9075 FORMAT(/10X,'THE MAX LAGRANGIAN FOR FMO PROJECTED ',
     *        'OUT ORBITALS IS ',F12.8/)
 9077 FORMAT(/10X,'THE MAX LAGRANGIAN FOR FMO PROJECTED ',
     *        'OUT ORBITALS WAS ',F12.8/)
 9080 FORMAT(1X,'EXCESSIVE NUMBER OF ITERATIONS...')
 9090 FORMAT(1X,'MCSCF IS NOT CONVERGED!'/
     *       1X,'A $VEC GROUP OF CURRENT MO-S IS IN THE PUNCH FILE'/
     *       1X,'USE THIS WITH GUESS=MOREAD TO RESTART THIS RUN')
 9100 FORMAT(/1X,'FINAL MCSCF ENERGY IS',F20.10,' AFTER',I4
     *       ,' ITERATIONS')
 9110 FORMAT(1X,'.....DONE WITH MCSCF ITERATIONS.....')
 9120 FORMAT(//10X,34(1H-)/10X,'CI BASED ON MCSCF NATURAL ORBITALS'/
     *         10X,34(1H-))
 9135 FORMAT(/10X,'INDUCED DIPOLE',5X,'ATOMIC UNITS',10X,'DEBYE',
     *       /10X,'      X(IND)  ',5X,F12.5,5X,F12.5,
     *       /10X,'      Y(IND)  ',5X,F12.5,5X,F12.5,
     *       /10X,'      Z(IND)  ',5X,F12.5,5X,F12.5)
C
 1010 FORMAT(1X,'CALCULATED ',F10.5,' THEORETICAL ',
     * F10.5,' NOT RENORMALIZED',F10.5)
 1100 FORMAT(1X,'CALCULATED ',F10.5,' THEORETICAL ',
     * F10.5,' RENORMALIZED',F10.5)
 1110 FORMAT(1X,'CALCULATED',F10.5,' (THEOR= ',
     * F10.5,')  ESCAPED',F10.5,' FINAL', F10.5)
      END
C*MODULE MCSCF   *DECK NEWTON
      SUBROUTINE NEWTON(CISTEP,NPFLG8,METHOD,ACURCY,DAMP,EKT,
     *                  CANONC,NPUNCH,FCORE,FORS,NOCI,
     *                  NWORD,NORB,NOROT,MXNORO,IROTMX,JROTMX,
     *                  MOFRZ,FINAL,MCNOCI,DOEXCH,DROPC,FULLNR)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION METHOD
C
      PARAMETER (MXATM=2000, MXRT=100, MXAO=8192)
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH)
C
      LOGICAL EKT,CANONC,CANSAV,FCORE,FORS,NOCI,FINAL,MCNOCI,
     *        DBGGS,DBUG,OUT,SOME,GOPARR,DSKWRK,MASWRK,SVDSKW,
     *        PACK2E,SOMSAV,MEM,NOTRF,DIRTRF,SCHWRZ,DOEXCH,DROPC,
     *        JACOBI,DUPAO,HSSEND,FULLNR,FDIRCT,QCORR
C
      DIMENSION NOROT(2,MXNORO),MOFRZ(15)
C
      COMMON /CIFILS/ NFT11,NFT12,NFT13,NFT14,NFT15,NFT16,IDAF20,NEMEMX
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,SDET,SZDET,
     *                GRPDET,STSYM,GLIST,NFLGDM(MXRT),
     *                IWTS(MXRT),NCORSV,NCORDET,NACTDET,NORBDET,
     *                NADET,NBDET,KDET,KSTDET,IROOT,IPURES,MAXW1,NITDET,
     *                MAXP,NCIDET,IGPDET,KSTSYM,NFTGCI
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FCCWFN/ NSPACE,MSTA(51),MNUM(51),MINI(51),MAXI(51),
     *                IAMI(51),IAMA(51),IBMI(51),IBMA(51),IDIM(51),
     *                LBST(51),NREF0,FDIRCT,QCORR,C0SQ
      COMMON /FMCOM / X(1)
      COMMON /HSSCTL/ IHESSM,IHREP,HSSEND
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /INT2IC/ NINTIC,ININTIC,NXXIC,LBUFPIC,LIXIC,LABSIX,NINTIX
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /JACOBI/ JACOBI,NJAOR,ELAST,ISTAT
      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 /NTNOPT/ ENERGY,ENERG0,DEMAX,SQCDF,ITER,MICRO,NOTRF
      COMMON /ORDOPT/ IDUM(6),NOSQUR
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
      COMMON /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FE,UNZOR,PORT
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IKREP,IEF,IP_F,NFMOPCM
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      INTEGER         D_OOOO,D_VOOO,D_VVOO,D_VOVO,D_VVVO,D_VVVV,
     *                D_OOOOAB,D_OOOOBB,D_VOOOAB,D_VOOOBA,D_VOOOBB,
     *                D_VVOOAB,D_VVOOBA,D_VVOOBB,D_VOVOAB,D_VOVOBB,
     *                D_U,D_UB,D_E,D_EB
      LOGICAL         NDOOOO,NDVOOO,NDVVOO,NDVOVO,NDVVVO,NDVVVV,NDCORE,
     *                NDVVOOBA,NDVVOOAB,NDVVOOBB,NDVOVOAB,NDVOVOBB,
     *                NDVOOOBA,NDVOOOAB,NDVOOOBB,NDOOOOAB,NDOOOOBB
      COMMON /TRFDMS/ D_OOOO,D_VOOO,D_VVOO,D_VOVO,D_VVVO,D_VVVV,
     *                D_OOOOAB,D_OOOOBB,D_VOOOAB,D_VOOOBA,D_VOOOBB,
     *                D_VVOOAB,D_VVOOBA,D_VVOOBB,D_VOVOAB,D_VOVOBB,
     *                D_U,D_UB,D_E,D_EB,
     *                NDOOOO,NDVOOO,NDVVOO,NDVOVO,NDVVVO,NDVVVV,NDCORE,
     *                NDVVOOBA,NDVVOOAB,NDVVOOBB,NDVOVOAB,NDVOVOBB,
     *                NDVOOOBA,NDVOOOAB,NDVOOOBB,NDOOOOAB,NDOOOOBB
C
      PARAMETER (TOL=5.0D-02, ZERO=0.0D+00, HALF=0.5D+00)
C
      DATA DM2,TEI/8HDM2     ,8HTEI     /
      DATA CHECK/8HCHECK   /, DEBUG/8HDEBUG   /, DBUGME/8HNEWTON  /
      DATA OPT,SAD/8HOPTIMIZE,8HSADPOINT/
      DATA RHF/8HRHF     /
      DATA GUGA,ALDET,GENCI,ORMAS,GMCCI
     *      /8HGUGA    ,8HALDET   ,8HGENCI   ,8HORMAS   ,8HGMCCI   /
      DATA RNONE/8HNONE    /
C
C     ----- DRIVER FOR NEWTON-RAPHSON ORBITAL IMPROVMENT -----
C     BUT THE BOTTOM HALF OF THE ROUTINE IS CALLED BY ALL CONVERGERS
C     IN ORDER TO DO THE FINAL PRINTOUT OF ORBITALS, CI VECTOR, ...
C
      DO 10 I=1,MXAO
         IA(I)=(I*I-I)/2
   10 CONTINUE
C
      DBGGS = (NPFLG8.GT.0 .OR. EXETYP.EQ.DEBUG .OR. EXETYP.EQ.DBUGME)
     *        .AND.  GOPARR
      OUT   = (NPFLG8.GT.0 .OR. EXETYP.EQ.DEBUG .OR. EXETYP.EQ.DBUGME)
     *        .AND. MASWRK
      DBUG  =   NPFLG8.GE.2
     *        .AND. MASWRK
      SOME  = ((NPFLG8.NE.-5.AND.ITER.EQ.1) .OR. OUT .OR. DBUG)
     *        .AND. MASWRK
      IF(SOME.AND..NOT.FINAL) WRITE(IW,9000)
C
C     ----- OBTAIN MAXIMUM SIZE OF FAST MEMORY -----
C     NOTE THAT THE MEMORY ALLOCATION HERE IS A SERIES OF GETFM
C     AND RETFM CALLS, BUT THAT MOST OF THE TIME RETFM IS CALLED,
C     THE CODE IS NOT DONE WITH ALL OF THE ARRAYS IT HAS BEEN USING.
C     THESE ARE REALLOCATED AT THE SAME POINTERS INTO THE MEMORY
C     BLOCK BEFORE THE NEXT GETFM, TO REUSE THESE ARRAYS.
C
      CALL GOTFM(NGOT)
      IF(NWORD.GT.0) NGOT=MIN(NWORD,NGOT)
C
C     ----- LEARN NUMBER OF CORE AND OCCUPIED ORBITALS -----
C
      IF(CISTEP.EQ.ALDET  .OR.
     *   CISTEP.EQ.GENCI  .OR.
     *   CISTEP.EQ.ORMAS  .OR.
     *   CISTEP.EQ.GMCCI) THEN
         NORBS = NCORSV+NACTDET
         NCORBS= NCORSV
      END IF
      IF(CISTEP.EQ.GUGA) THEN
         SVDSKW=DSKWRK
         DSKWRK=.TRUE.
         CALL SEQREW(NFT11)
         CALL RDDRT(NORBMX,NORBS,NSYM,NROWS,NWKS,LEVFRM,NEMEMX,
     *              NREFS,IEXCT,NFOCI,INTACT,NCORBS)
         CALL SEQREW(NFT11)
         DSKWRK=SVDSKW
         IF (DROPC) NORBS = NORBS + NCORBS
      END IF
C
C     ----- AT THIS POINT IT IS GOOD TO REMEMBER THAT
C     NORB   = TOTAL NUMBER OF MOLECULAR ORBITALS (OCC + VIRT),
C              BEING USED IN THE MCSCF (NORMALLY NUM).
C     NORBS  = NUMBER OF OCCUPIED MOLECULAR ORBITALS (CORE + ACT)
C     NCORBS = NUMBER OF CORE ORBITALS (CLOSED SHELL)
C
      M1=NORB
      M2=(M1*(M1+1))/2
      M4=(M2*(M2+1))/2
      N0=NQMT
      N1=NUM
      N2=(N1*(N1+1))/2
      N3=N1*N1
      L1=NORBS
      L2=(L1*(L1+1))/2
      NLGR=NORB*NORBS
      NCI=NORBS-NCORBS
      LG1=(NCI*(NCI+1))/2
C
C     ----- FINAL PASS FOR PRINTOUT -----
C
      IF(FINAL) GO TO 500
C
C     ----- SET UP ORBITAL PAIRS ROTATION INDICES -----
C     NOTE THAT MEMORY IS NOT FORMALLY ALLOCATED AT THIS POINT.
C
      CALL VALFM(LOADFM)
      I00 = LOADFM + 1
      I10 = I00 + NLGR
      I11 = I10 + NLGR
      I12 = I11 + L1
C
      CALL NTNPAR(NROT,NROT2,X(I00),X(I10),X(I11),NOROT,MXNORO,MOFRZ,
     *            NORB,NORBS,NCORBS,FORS,FCORE,DBUG)
C
C  DISTRIBUTED DATA FULL NEWTON-RAPHSON MCSCF WITHOUT
C  STORING THE AUGMENTED HESSIAN MATRIX
C  WORKS FOR DROPC=.TRUE. ONLY
C
      IF (GOPARR) THEN
        IF(SOME) WRITE(IW,9900)
        NBF  = NUM
        NMOS = NORB
        NOCC = NORBS
        NCOR = NCORBS
        NACT = NOCC - NCOR
        NVIR = NMOS - NOCC
        MXPN = 50
C
        NATR = (NACT*NACT+NACT)/2
        NMTR = (NMOS*NMOS+NMOS)/2
        NAT4 = (NATR*NATR+NATR)/2
        MXTR = (MXPN*MXPN+MXPN)/2
        NVSQ = NVIR*NVIR
        MXSQ = MXPN*MXPN
C
C  DATA STRUCTURE LENGTHS
C
        LIHES = NMOS*NOCC
        LIROT = NMOS*NOCC
        LLAGN = NMOS*NOCC
        LOPDM = NATR
        LTPDM = NAT4
        LFCOR = NMTR
        LFVAL = NMTR
        LDIAH = NROT+1
        LBUFF = MAX(NVSQ,2*NINTMX)
        LTVEC = MXPN*(NROT+1)
        LPROD = MXPN*(NROT+1)
        LCVEC = NROT+1
        LDMTX = MXTR
        LSDMX = MXTR
        LSDVC = MXSQ
        LSDEI = MXPN
        LSDWS = MXPN*8
        LSDIN = MXPN
C
C  COMPUTE ADDRESSES, VALFM CALLED ONCE ABOVE.
C
        IIROT = LOADFM + 1       ! INDEX OF ROTATIONS TO OPTIMIZE
        IIHES = IIROT  + LIROT   ! INDEX OF NON-REDUNDANT ROTATIONS
        ILAGN = IIHES  + LIHES   ! LAGRANGIAN
        IOPDM = ILAGN  + LLAGN   ! 1-EL DENSITY OVER ACTIVE MOS
        ITPDM = IOPDM  + LOPDM   ! 2-EL DENSITY OVER ACTIVE MOS
        IFCOR = ITPDM  + LTPDM   ! CORE FOCK MATRIX
        IFVAL = IFCOR  + LFCOR   ! VALENCE FOCK MATRIX
        IDIAH = IFVAL  + LFVAL   ! DIAGONAL OF AUGMENTED HESSIAN
        IBUFF = IDIAH  + LDIAH   ! MESSAGE, I/O BUFFER
        ITVEC = IBUFF  + LBUFF   ! TRIAL VECTORS
        IPROD = ITVEC  + LTVEC   ! TRIAL VECTOR PRODUCTS
        ICVEC = IPROD  + LPROD   ! CORRECTION VECTOR
        IDMTX = ICVEC  + LCVEC   ! DAVIDSON MATRIX
        ISDMX = IDMTX  + LDMTX   ! SMALL DIAG MATRIX
        ISDVC = ISDMX  + LSDMX   ! SMALL DIAG VECTORS
        ISDEI = ISDVC  + LSDVC   ! SMALL DIAG EIGENVALUES
        ISDWS = ISDEI  + LSDEI   ! SMALL DIAG WORKSPACE
        ISDIN = ISDWS  + LSDWS   ! SMALL DIAG INDEX
        LAST  = ISDIN  + LSDIN
        NEED  = LAST   - LOADFM
        CALL GETFM(NEED)
        MXNEED = NEED
        IF (EXETYP.NE.CHECK) THEN
C
C  READ 1-PARTICLE DENSITY OVER ACTIVE MOS
C
          CALL DAREAD(IDAF,IODA,X(IOPDM),LOPDM,320,0)
C
C  READ 1-EL (CORE-HAMILTONIAN) INTEGRALS OVER MOS
C
          IHAMO = IFCOR
          CALL DAREAD(IDAF,IODA,X(IHAMO),NMTR,355,0)
C
C  COMPUTE 1-EL (CORE-HAMILTONIAN) ENERGY
C
          CALL CHNRGY(NCOR,NACT,X(IHAMO),X(IOPDM),E1EL)
          EC = E1EL
C
C  FORM VALENCE FOCK OPERATOR IN PARALLEL
C
          CALL FVADDI(NMOS,NCOR,NACT,X(IOPDM),X(IFVAL),X(IBUFF))
C
C  READ CORE FOCK MATRIX IN MO BASIS
C
          CALL DAREAD(IDAF,IODA,X(IFCOR),NMTR,35,0)
C
C  READ 2-EL DENSITY OVER ACTIVE MOS ON MASTER
C
          IINTS = IBUFF
          ILABS = IBUFF + NINTMX
          CALL DCOPY(NAT4,ZERO,0,X(ITPDM),1)
          CALL RD2PDM(NFT15,X(ITPDM),X(IINTS),X(ILABS),NINTMX)
C
C  BROADCAST 2-PDM TO ALL NODES
C
          CALL DDI_BCAST(2007,'F',X(ITPDM),NAT4,MASTER)
C
C  FORM LAGRANGIAN, TOTAL ENERGY
C
          CALL LAGDDI(NMOS,NCOR,NOCC
     *,               X(IOPDM),X(ITPDM),X(IFCOR),X(IFVAL)
     *,               X(IBUFF),X(ILAGN),ENERGY)
          ENERGY = ENERGY + EC*HALF + ENUCR
C
C  SAVE LAGRANGIAN (AS AN N**2 ARRAY?)
C
          CALL DAWRIT(IDAF,IODA,X(ILAGN),NBF*NBF,403,0)
C
C  GET LAGRANGIAN ASYMMETRY, USE IT TO DETERMINE WHICH
C  ROTATIONS NEED OPTIMIZING...CREATE A NEW INDEX FOR
C  THESE (IROT), ADJUST NUMBER OF ROTATIONS (NROT).
C
          NSKIP = 0
          IF (FCORE) NSKIP = NCOR
          CALL NTNROT(DEMAX,IROTMX,JROTMX,X(ILAGN),X(IIROT)
     *,               X(IIHES),NROT,NORB,NORBS,NSKIP,ACURCY,SOME,OUT)
C
C  WRITE ORBITAL ROTATION TABLE ON DAF RECORD 322
C
          CALL DAWRIT(IDAF,IODA,X(IIROT),LIROT,322,1)
C
C  DIAGONAL ELEMENTS OF THE AUGMENTED HESSIAN
C
          CALL DIADDI(NMOS,NCOR,NOCC,NROT
     *,               X(IOPDM),X(ITPDM),X(IFCOR),X(IFVAL)
     *,               X(ILAGN),X(IIROT),X(IBUFF),X(IDIAH))
C
C  FIND THE LOWEST ROOT OF THE AUGMENTED HESSIAN USING
C  DAVIDSON METHOD. CONTRIBUTIONS TO THE HESSIAN ARE MULTIPLIED
C  BY ELEMENTS OF THE TRIAL VECTOR AND SUMMED DIRECTLY TO THE
C  PRODUCT-VECTOR. IN THE PARALLEL CALL TO NTNDVD, THE STORAGE
C  ADDRESS OF THE AUGMENTED HESSIAN IS SET TO 'DUMMY'.
C
          MXIT  = 200
          CGTL  = ACURCY
          PRTL  = 1.0D-09   ! A DEBUG SETTING
          DBUG  = .FALSE.
          CALL NTNDVD(DUMMY
     *,               X(IDIAH),X(ITVEC),X(IPROD),X(ICVEC),X(IDMTX)
     *,               X(ISDMX),X(ISDVC),X(ISDEI),X(ISDWS),X(ISDIN)
     *,               IA,NROT+1,MXPN,CGTL,MXIT,OUT,DBUG,PRTL
     *,               NMOS,NCOR,NOCC,NROT  ! FOR COMPUTING AUG HESSIAN
     *,               X(IOPDM),X(ITPDM),X(IFCOR),X(IFVAL)
     *,               X(ILAGN),X(IIROT),X(IBUFF))
          IF (SOME) WRITE(IW,9560)
          IF (SOME) CALL TIMIT(1)
C
C  LENGSFIELD'S EMPIRICAL FUDGING
C
          ITEMP = IPROD
          CALL NTNLEF(X(ITVEC),X(ITEMP),X(ILAGN),X(IIROT)
     *,               NROT,NMOS,NOCC,OUT,SQCDF)
C
C  REMOVE DISTRIBUTED MATRICES IN REVERSE ORDER TO CREATE'S
C
          CALL DDI_DESTROY(D_VVOO)
          CALL DDI_DESTROY(D_OOOO)
          CALL DDI_DESTROY(D_VOOO)
          CALL DDI_DESTROY(D_VOVO)
        END IF  ! EXETYP
C
C  NOTE: IT IS CUSTOMARY (HERE) TO FORMALLY FREE THE MEMORY
C        THEN RE-USE CERTAIN ADDRESSES FOR THE NEXT STEP.
C
        CALL RETFM(NEED)
      ELSE
C
C   END OF PARALLEL SECTION, AND TOP OF SEQUENTIAL SECTION.
C
      IF (.NOT. DROPC) GO TO 40
C
C     ----- FORM VALENCE FOCK OPERATOR -----
C
      NIA=MAX0(NROT,N2)
      I12= I11+NLGR
      I15= I12+NROT2
      I20= I15+LG1
      I30= I20+N2
      I40= I30+N2
      I50= I40+N2
      I55= I50 + N2
      I60= I55 + N1
      I70= I60+MAX0(NINTMX,N3)
      I80= I70+MAX0(NINTMX,N2)
      LAST= I80 + NIA
      I60A=I60
      I70A=I70
      IF(NINTIC.NE.0) THEN
         I60A=LBUFPIC
         I70A=LIXIC
      ENDIF
      IF(DIRTRF) THEN
         NSH2 = (NSHELL*NSHELL+NSHELL)/2
         CALL BASCHK(LMAX)
                       NANGM =  4
         IF(LMAX.EQ.2) NANGM =  6
         IF(LMAX.EQ.3) NANGM = 10
         IF(LMAX.EQ.4) NANGM = 15
         IF(LMAX.EQ.5) NANGM = 21
         IF(LMAX.EQ.6) NANGM = 28
         MAXG = NANGM**4
         LXINTS= LAST
         LGHOND= LXINTS + NSH2
         LDSH  = LGHOND + MAXG
         LDDIJ = LDSH   + NSH2
         LAST  = LDDIJ  + 49*MXG2
      ELSE
         LXINTS= LAST
         LGHOND= LAST
         LDSH  = LAST
         LDDIJ = LAST
      END IF
      NEED=LAST-LOADFM-1
      IF(SOME  .AND.  NEED.GT.NGOT) WRITE(IW,9730) NEED,NGOT
      CALL GETFM(NEED)
C
      IF(EXETYP.EQ.CHECK) GO TO 30
      IF(SOME) WRITE(IW,9740) NEED
C
C     ----- READ 1E DENSITY -----
C
      CALL DAREAD(IDAF,IODA,X(I15),LG1,320,0)
C
      CALL NTNDE1(X(I60),X(I20),X(I30),X(I15),X(I80),NORBS,NCORBS,
     *            X(I50),EC,NIA)
C
      IF(DIRTRF) THEN
         NINT=0
         NSCHWZ=0
         SCHWRZ = ISCHWZ.GT.0
         DUMMY = ZERO
         IF(DOEXCH) THEN
            CALL EXCHNG(X(LXINTS),X(LGHOND),X(LDDIJ),NSH2,MAXG,INTTYP)
         ELSE
            CALL DAREAD(IDAF,IODA,X(LXINTS),NSH2,54,0)
         END IF
         IF(SCHWRZ) CALL SHLDEN(RHF,X(I30),DUMMY,DUMMY,X(LDSH),
     *                                         IA,N1,N2,NSH2,1)
         CALL VCLR(X(I50),1,N2)
         CALL TWOEI(RHF,.TRUE.,.FALSE.,.FALSE.,.FALSE.,
     *              INTTYP,SCHWRZ,NINT,NSCHWZ,N1,N2,
     *              DUMMY,DUMMY,IDUMMY,1,
     *              X(LXINTS),NSH2,X(LGHOND),MAXG,X(LDDIJ),
     *              IA,X(I30),X(I50),DUMMY,DUMMY,X(LDSH),
     *              DUMMY,DUMMY,1)
         CALL DSCAL(N2,HALF,X(I50),1)
         II=I50-1
         DO 110 I=1,N1
            II = II+I
            X(II) = X(II) + X(II)
  110    CONTINUE
      ELSE
         SVDSKW = DSKWRK
         DSKWRK = .TRUE.
         CALL SEQREW(IJKO)
         DUPAO  = ITRFAO.EQ.1
         CALL HSTAR(X(I30),X(I50),X(I60A),X(I70A),NINTMX,IA,NOPK,DUPAO)
         DSKWRK = SVDSKW
      END IF
C
C     ----- ADD PARTIAL FOCK MATRICES IF DISTRIBUTED PROCESSING
C
      IF (GOPARR) CALL DDI_GSUMF(2200,X(I50),N2)
C
C     ----- SYMMETRIZE 'SKELETON' FOCK MATRICES -----
C
      CALL SYMH(X(I50),X(I70),IA)
C
C     ----- TRANSFORMATION TO MO'S -----
C
      CALL DAREAD(IDAF,IODA,X(I40),N2,35,0)
      CALL DAREAD(IDAF,IODA,X(I60),N3,15,0)
C
      CALL VCLR(X(I55),1,N1)
      CALL TFTRI(X(I20),X(I40),X(I60),X(I55),NORB,NUM,NUM)
      CALL TFTRI(X(I30),X(I50),X(I60),X(I55),NORB,NUM,NUM)
C
C     ----- FOCK CONTRIBUTIONS TO LAGRANGIAN AND HESSIAN -----
C
C           ----- POINTERS
C
C       ADRESS  NAME  LENGTH COMMENT
C
C       AT I00   IB    NLGR    INDEX FOR THE LAGRANGIAN
C       AT I10   IC    NLGR    INDEX FOR THE HESSIAN
C       AT I11   EG    NLGR    LAGRANGIAN MATRIX
C       AT I12   EH    NROT2   HESSIAN MATRIX
C       AT I15  GIJ    LG1     DM1
C       AT I20  FFC    N2      CORE FOCK MATRIX OVER MO'S
C       AT I30  FFA    N2      VALENCE FOCK MATRIX OVER MO'S
C       AT I80  IA     NIA     ROTATION INDEX
C     NOTE POINTERS I00, I10, I11, I12, I15 ARE REPEATED BELOW IN
C     THE COMPUTATIONS INVLOVING TRANSFORMED INTEGRALS.
C
      CALL VCLR(X(I11),1,NLGR)
      CALL VCLR(X(I12),1,NROT2)
C
      IF(MASWRK) THEN
         CALL NTNLAG(X(I80),X(I20),X(I30),X(I15),X(I11),
     *               NORB,NORBS,NCORBS,X(I00),NIA)
         CALL NTNHES(X(I80),X(I10),X(I20),X(I30),X(I15),X(I12),
     *               NORB,NORBS,NCORBS,NIA)
      END IF
C
      IF(SOME) CALL TIMIT(1)
   30 CONTINUE
      CALL RETFM(NEED)
C
C     ----- LEARN MEMORY DEMANDS FOR THE 3 METHODS -----
C     DON'T WRITE THIS NEXT CODE AS AN IF-THEN, SOME COMPILERS
C     DON'T LIKE BRANCHING BACK INTO AN IF-THEN BLOCK.
C
   40 CONTINUE
      MEM = .FALSE.
      SOMSAV=SOME
      MXNEED=0
      IF(DROPC) MXNEED = NEED
      NEED1=0
      NEED2=0
      NWKS=0
      IF((ITER.GT.1) .OR. (NEVALS.GT.0)) GO TO 90
         MEM=.TRUE.
         SOME=.FALSE.
         MXNEED=0
         GO TO 100
   50    NEED1=MXNEED
         MXNEED=0
         GO TO 200
   60    NEED2=MXNEED
         MXNEED=0
         MEM=.FALSE.
         SOME=SOMSAV
         IF (MASWRK) WRITE(IW,9010) NEED1,NEED2,NGOT
         IF(METHOD.EQ.DM2) NEED = NEED1
         IF(METHOD.EQ.TEI) NEED = NEED2
         IF(MASWRK  .AND.  NEED.GT.NGOT)
     *      WRITE(IW,*) 'INSUFFICIENT MEMORY, 2E- TERMS NEED=',NEED
   90 CONTINUE
C
C     ----- CHOOSE METHOD FOR 2 ELECTRON NTN-RPH TERMS -----
C
      IF(EXETYP.EQ.CHECK) THEN
         CALL GETFM(NEED)
         CALL RETFM(NEED)
         GO TO 400
      END IF
      IF(METHOD.EQ.DM2)    GO TO 100
      IF(METHOD.EQ.TEI)    GO TO 200
      IF(MASWRK) WRITE(IW,*) 'UNKNOWN HESSIAN METHOD'
      CALL ABRT
C
C           ****************
C     ***** THIS IS METHOD 1 *****
C           ****************
C
C     CALCULATE ORBITAL HESSIAN DIRECTLY WITHOUT FORMULAS.
C     THIS METHOD IS DRIVEN BY THE SECOND ORDER DENSITY MATRIX.
C
  100 CONTINUE
      IF(SOME) WRITE(IW,9100)
C
C     ----- FORM 2E-CONTRIBUTION TO LAGRANGIAN + HESSIAN MATRIX -----
C
      MI    = NORBS
      MIJ   =(MI*(MI+1))/2
      MIJKL =(MIJ*(MIJ+1))/2
      MA    =(NORB-NORBS)
      MAJ   =MA*MI
      MAJKL =MAJ*MIJ
      MAB   =(MA*(MA+1))/2
      MABKL =MAB*MIJ
      MAJBL =(MAJ*(MAJ+1))/2
      MAXINT= MAX(MIJKL,MAJKL,MABKL,MAJBL)
      NIA=(NORB*(NORB+1))/2
C
C     ORB. ROT. INDEX     AT   IX(I10)
C     LAGRANGIAN MATRIX   AT    X(I11)
C     HESSIAN MATRIX      AT    X(I12)
C     DM1                 AT    X(I15)
C     ORDERED -INT2-      AT    X(I20)
C     -INT2- BUFFER       AT    X(I30)
C     -INT2- LABELS       AT   IX(I40)
C     PAIR INDEX          AT   IX(I50)
C
C     NOTE POINTERS I00, I10, I11, I12, I15 ARE INHERITED FROM ABOVE
C
      I00=LOADFM+1
      I10=I00+NLGR
      I11=I10+NLGR
      I12=I11+NLGR
      I15=I12+NROT2
      I20=I15+LG1
      I30=I20+MAXINT
      I40=I30+NINTMX
      I50=I40+NINTMX*LABSIZ
      I60=I50+NIA
      NEED=I60-I00
      MXNEED=MAX(NEED,MXNEED)
      IF(MEM) GO TO 50
      IF(NEED.GT.NGOT) THEN
         IF (MASWRK) WRITE(IW,9110) NEED,NGOT
         CALL ABRT
      END IF
      CALL GETFM(NEED)
C
      NFTG=NFT15
      NFTX=IJKT
C
C     ----- -LGR- + -COU- + -EXC- FROM (IJ//KL) -----
C
      IF(SOME) WRITE(IW,9120)
      CALL SEQREW(NFTX)
      IF (MASWRK) READ(NFTX)
      CALL M1EGHR(NFTX,X(I20),X(I30),X(I40),X(I50),NINTMX,MI,NOSQUR)
      CALL M1EGH(NLGR,NROT2,X(I11),X(I12),X(I20),NFTG,X(I30),X(I15),
     *         X(I40),NIA,X(I50),NORB,X(I00),X(I10),NINTMX,MI,NCORBS,
     *         DROPC)
      IF(SOME) CALL TIMIT(1)
      IF(MA.EQ.0) CALL RETFM(NEED)
      IF(MA.EQ.0) GO TO 400
C
C     ----- -LGR- + -COU- + -EXC- FROM (AJ//KL) -----
C
      IF(SOME) WRITE(IW,9130)
      CALL SEQREW(NFTX)
      IF (MASWRK) READ(NFTX)
      CALL M1LGRR(NFTX,X(I20),X(I30),X(I40),X(I50),
     *            NINTMX,MI,MA,MIJ,NOSQUR)
      CALL M1LGR(NLGR,NROT2,X(I11),X(I12),X(I20),NFTG,X(I30),X(I15),
     *           X(I40),NIA,X(I50),NORB,X(I00),X(I10),
     *           NINTMX,MI,MA,MIJ,NCORBS,DROPC)
      IF(SOME) CALL TIMIT(1)
C
C     ----- -COU- FROM (AB//KL) -----
C
      IF(SOME) WRITE(IW,9140)
      CALL SEQREW(NFTX)
      IF (MASWRK) READ(NFTX)
      CALL M1COUR(NFTX,X(I20),X(I30),X(I40),X(I50),
     *            NINTMX,MI,MA,MIJ,MAB,NOSQUR)
      CALL M1COU(NROT2,X(I12),X(I20),NFTG,X(I30),X(I15),X(I40),
     *           NIA,X(I50),NORB,X(I10),NINTMX,MI,MA,MAB,NCORBS,
     *           DROPC)
      IF(SOME) CALL TIMIT(1)
C
C     ----- -EXC- FROM (AJ//BL) -----
C
      IF(SOME) WRITE(IW,9150)
      CALL SEQREW(NFTX)
      IF (MASWRK) READ(NFTX)
      CALL M1EXCR(NFTX,X(I20),X(I30),X(I40),X(I50),
     *            NINTMX,MI,MAJ,NOSQUR)
      CALL M1EXC(NROT2,X(I12),X(I20),NFTG,X(I30),X(I15),X(I40),
     *           NIA,X(I50),NORB,X(I10),NINTMX,MI,MA,NCORBS,DROPC)
C
      CALL RETFM(NEED)
      GO TO 400
C
C           ****************
C     ***** THIS IS METHOD 2 *****
C           ****************
C
C     CALCULATE ORBITAL HESSIAN DIRECTLY WITHOUT FORMULAE.
C     METHOD DRIVEN BY 2-ELECTRON INTEGRALS.
C           THIS METHOD DOES NOT REQUIRE MUCH CORE BEYOND
C           WHAT IS NEED FOR THE -GRADIENT- + HESSIAN, BUT
C           IT IS ...... SLOWER THAN METHODS 1 AND 3.
C
  200 CONTINUE
      IF(SOME) WRITE(IW,9200)
C
C     ----- FORM 2E-CONTRIBUTION TO LAGRANGIAN + HESSIAN MATRIX -----
C
      M2=(NORBS*(NORBS+1))/2
      M4=(M2*(M2+1))/2
      NIA=(NORB*(NORB+1))/2
C
C     NOTE POINTERS I00, I10, I11, I12 ARE INHERITED FROM ABOVE
C
      I00=LOADFM+1
      I10=I00+NLGR
      I11=I10+NLGR
      I12=I11+NLGR
      I20=I12+NROT2
      I30=I20+M4
      I40=I30+NINTMX
      I50=I40+NINTMX*LABSIZ
      I60=I50+NIA
      NEED=I60-I00
      MXNEED=MAX(NEED,MXNEED)
      IF(MEM) GO TO 60
      IF(NEED.GT.NGOT) THEN
         IF (MASWRK) WRITE(IW,9210) NEED,NGOT
         CALL ABRT
      END IF
      CALL GETFM(NEED)
C
C     ORB. ROT. INDEX     AT   IX(I10)
C     LAGRANGIAN MATRIX   AT    X(I11)
C     HESSIAN MATRIX      AT    X(I12)
C     ORDERED -INT2-      AT    X(I20)
C     -INT2- BUFFER       AT    X(I30)
C     -INT2- LABELS       AT   IX(I40)
C     PAIR INDEX          AT   IX(I50)
C
      NFTG=NFT15
      CALL M2DM2R(NFTG,NORBS,X(I20),X(I30),X(I40),NINTMX,X(I50))
      NFTX=IJKT
      CALL SEQREW(NFTX)
      IF(MASWRK) READ(NFTX)
      CALL M2TEI(NLGR,NROT2,X(I11),X(I20),NFTX,X(I30),X(I40),
     *           X(I50),X(I00),X(I10),NINTMX,NIA,NORB,NORBS)
      CALL RETFM(NEED)
C
C     ***** ALL METHODS *****
C     ----- FORM 1E-CONTRIBUTION TO LAGRANGIAN AND HESSIAN -----
C
  400 CONTINUE
      IF(SOME) WRITE(IW,9500)
      IF(SOME) CALL TIMIT(1)
      NIA=MAX(NROT,NORB)
      M2=(NORB *(NORB +1))/2
      L2=(NORBS*(NORBS+1))/2
C
C     NOTE POINTERS I00, I10, I11, I12 ARE INHERITED FROM ABOVE
C
      I00=LOADFM+1
      I10=I00+NLGR
      I11=I10+NLGR
      I12=I11+NLGR
      I20=I12+NROT2
      I30=I20+M2
      I40=I30+L2
      I50=I40+NIA
      I60=I50+NORB
      NEED=I60-I00
      IF(NEED.GT.NGOT) THEN
         IF (MASWRK) WRITE(IW,9510) NGOT,NEED
         CALL ABRT
      END IF
      MXNEED=MAX(NEED,MXNEED)
      CALL GETFM(NEED)
C
C         DEBUG PRINTOUT OF MATRICES WITH 2E- TERMS COMPLETED ABOVE
C
      IF(DBGGS) THEN
         CALL DDI_GSUMF(2000,X(I11),NLGR)
         CALL GSUMCP(2010,X(I12),NROT2)
      END IF
      IF(OUT) THEN
         WRITE(IW,*) '2 ELECTRON PART OF LAGRANGIAN'
         CALL PRSQ(X(I11),NORB,NORBS,NORBS)
         WRITE(IW,*) '2 ELECTRON PART OF HESSIAN'
         CALL PRTRI(X(I12),NROT)
      END IF
      IF(DBGGS) THEN
         CALL DSCAL(NLGR ,1.0D+00/NPROC,X(I11),1)
         CALL DSCAL(NROT2,1.0D+00/NPROC,X(I12),1)
      END IF
C
C           ONE ELECTRON TERMS
C
      IF(EXETYP.EQ.CHECK) THEN
         CALL VCLR(X(I11),1,N3)
         ENERGY=ZERO
      ELSE
         IF(GOPARR  .OR.  DROPC) THEN
            CALL NTNRGY(X(I11),X(I12),X(I40),X(I50),
     *                  NLGR,NROT,NROT2,NORB,NORBS,
     *                  NIA,NFT16,ENERGY)
            ENERGY=ENERGY + (EC * HALF) + ENUCR
         ELSE
            CALL NTNONE(X(I11),X(I12),X(I20),X(I30),X(I40),X(I50),
     *                  X(I10),NLGR,NROT,NROT2,NORB,M2,NORBS,L2,
     *                  NIA,IJKT,NFT15,NFT16,ENERGY)
            ENERGY=ENERGY+ENUCR
         END IF
      END IF
C
      IF(OUT) THEN
         WRITE(IW,*) 'LAGRANGIAN MATRIX'
         CALL PRSQ(X(I11),NORB,NORBS,NORBS)
         WRITE(IW,*) 'HESSIAN MATRIX (1 AND 2 E- INTEGRAL TERMS)'
         CALL PRTRI(X(I12),NROT)
      END IF
      IF(SOME) WRITE(IW,9520)
      IF(SOME) CALL TIMIT(1)
C
C     ----- SAVE LAGRANGIAN, CHECK ITS ASYMMETRY -----
C
      CALL DAWRIT(IDAF,IODA,X(I11),N3,403,0)
C
C     IGNORE ANY CORE ROTATION, INCLUDING CORE-CORE FOR FCORE=.T..
C
      NSKIP=0
      IF(FCORE) NSKIP=NCORBS
C
      CALL NTNROT(DEMAX,IROTMX,JROTMX,X(I11),X(I00),X(I10),
     *            MROT,NORB,NORBS,NSKIP,ACURCY,SOME,OUT)
C
C     ----- WRITE ORBITAL ROTATION TABLE ON DAF RECORD 322 -----
C
      CALL DAWRIT(IDAF,IODA,X(I00),NLGR,322,1)
C
      IF(OUT) WRITE(IW,9530) ITER,MICRO,ENERGY,DEMAX
      SQCDF=ZERO
      IF(MROT.EQ.0) THEN
         CALL RETFM(NEED)
         GO TO 450
      END IF
      IF(MROT.GT.NROT) THEN
         IF (MASWRK) WRITE(IW,9540) MROT,NROT
         CALL ABRT
      END IF
      CALL RETFM(NEED)
C
      IF (NOTRF) GO TO 430
C
C     ----- SET UP AUGMENTED HESSIAN MATRIX AND GET LOWEST ROOT -----
C     POINTERS I00, I10, I11, I20 ARE INHERITED FROM ABOVE
C
      MROT1  = MROT+1
      MROT2  =(MROT1*MROT1+MROT1)/2
      MXXPAN = 50
      MXTR   = (MXXPAN*MXXPAN+MXXPAN)/2
      NIA    = MAX(MROT1,MXXPAN)
C
      I00  = LOADFM + 1                 !  LAGRANGIAN INDEX
      I10  = I00    + NLGR              !  ROTATION INDEX
      I11  = I10    + NLGR              !  LAGRANGIAN
      I12  = I11    + NLGR              !  HESSIAN
      I20  = I12    + MROT2             !  PAIR INDEX
      I30  = I20    + NIA               !  I/O BUFFER
      I40  = I30    + MROT              !  TEMPORARY
      I50  = I40    + MROT              !  HESSIAN DIAGONAL
      I51  = I50    + MROT1             !  DAVIDSON TRIAL VECTORS
      I52  = I51    + MROT1*MXXPAN      !  TRIAL VECTOR WORKSPACE
      I53  = I52    + MROT1*MXXPAN      !  WORKSPACE
      I70  = I53    + MROT1             !  DVDSN MATRIX
      I71  = I70    + MXTR              !  DVDSN WORKSPACE
      I72  = I71    + MXTR              !  DVDSN VECTORS
      I73  = I72    + MXXPAN*MXXPAN     !  DVDSN EIGENVALUES
      I74  = I73    + MXXPAN            !  DVDSN WORKSPACE
      I75  = I74    + MXXPAN*8          !  DVDSN INDEX
      LAST = I75    + MXXPAN
      NEED = LAST   - I00
      IF (NEED.GT.NGOT) THEN
         IF (MASWRK) WRITE(IW,9550) NGOT,NEED
         CALL ABRT
      END IF
      CALL GETFM(NEED)
      MXNEED = MAX(NEED,MXNEED)
C
      CALL NTNRPH(X(I00),X(I10),X(I11),X(I12),X(I20),X(I30),X(I40)
     *,           X(I50),X(I51),X(I52),X(I53)
     *,           X(I70),X(I71),X(I72),X(I73),X(I74),X(I75)
     *,           MROT,MROT1,MROT2,NIA,NORB,NORBS,NFT16,OUT
     *,           SQCDF,ACURCY,MXXPAN)
C
      CALL RETFM(NEED)
      IF(SOME) WRITE(IW,9560)
      IF(SOME) CALL TIMIT(1)
C
C          CLOSE PARALLEL/SEQUENTIAL IF/ELSE CLAUSES
C
      END IF
C
C          ARE WE CONVERGED?
C
      IF(DEMAX.LT.ACURCY) GO TO 450
C
C     ----- BUILD TRANSFORMATION MATRIX AND TRANSFORM ORBITALS -----
C
      I00=LOADFM+1
      I10=I00+NLGR
      I11=I10+NLGR
      I20=I11+NLGR
      I30=I20+NUM*NUM
      I40=I30+NORB*NORB
      I50=I40+NORB
      NEED=I50-I00
      IF(NEED.GT.NGOT) THEN
         IF (MASWRK) WRITE(IW,9570) NGOT,NEED
         CALL ABRT
      END IF
      MXNEED=MAX(NEED,MXNEED)
      CALL GETFM(NEED)
C
      CALL NTNTRF(NUM,NORB,NORBS,NCORBS,FCORE,DAMP,
     *            X(I11),X(I20),X(I30),X(I40))
      IF(OUT) THEN
         WRITE(IW,*) 'TRANSPOSED NEW ORBITALS IN OLD ORBITAL BASIS'
         CALL NTNOUT(X(I30),NORB,NORB)
         WRITE(IW,*) 'TRANSPOSED NEW ORBITALS IN -AO- BASIS'
         CALL NTNOUT(X(I20),NORB,NUM)
      END IF
      CALL RETFM(NEED)
      IF (.NOT. NOTRF) GO TO 450
C
C     ---- SET UP HESSIAN MATRIX AND SAVE ON FILE 23 FOR QUD-MCSCF ----
C
  430 CONTINUE
      SQCDF=ZERO
      IF (MROT.EQ.0) GO TO 450
      MROT1= MROT
      MROT2=(MROT1*(MROT1+1))/2
      NIA=MROT1
C
      I00=LOADFM+1
      I10=I00+NLGR
      I11=I10+NLGR
      I12=I11+NLGR
C
      I20=I12
      I30=I20+MROT2
      I40=I30+MROT1
      I50=I40+NIA
      I60=I50+NROT
      NEED=I60-I00
      IF(NEED.GT.NGOT) THEN
         IF (MASWRK) WRITE(IW,9550) NGOT,NEED
         CALL ABRT
      END IF
      CALL GETFM(NEED)
      CALL QUDHSS(X(I00),X(I10),X(I11),X(I20),X(I30),X(I40),X(I50),
     *            MROT,MROT1,MROT2,NIA,NORB,NORBS,NFT16,OUT)
C
C     ---- WRITE ORBITAL GRADIENT AND HESSIAN ON SEQUENTIAL FILE 23 ----
C
      CALL SEQREW(23)
      CALL SQWRIT(23,X(I30),MROT1)
      CALL SQWRIT(23,X(I20),MROT2)
C
      CALL RETFM(NEED)
C
C     ----- DONE WITH NTN-RPH ORBITAL IMPROVEMENT -----
C
  450 CONTINUE
      IF(SOME) WRITE(IW,9580)
      IF(SOME) CALL TIMIT(1)
      IF(SOME) WRITE(IW,9590) MXNEED
      RETURN
C
C     ----- FINAL=.TRUE. -----
C     PRINT WAVEFUNCTION, SET UP FOR PROPERTIES, ETC.
C
  500 CONTINUE
      IF(EXETYP.EQ.CHECK) RETURN
      SOME=NPRINT.NE.-5 .AND. MASWRK
      CALL VALFM(LOADFM)
C
C          I SUSPECT THIS SHOULD BE IN THE MCSCF DRIVER WITH ALL
C          THE OTHER FINAL SOLVATION ENERGY CORRECTIONS/PRINTING.
C
      IF(NFMOPCM.NE.0.AND.(CISTEP.EQ.ALDET  .OR.
     *                     CISTEP.EQ.GENCI  .OR.
     *                     CISTEP.EQ.GMCCI)) THEN
         DO II=1,KDET
            ESTATE(II)=ESTATE(II)+0.5D+00*(PC+PX)
         ENDDO
      END IF
C
C     ----- PRINT CI EXPANSION COEFFICIENTS -----
C
      IF(SOME) THEN
         IF(MCNOCI) THEN
            WRITE(IW,9700)
         ELSE
            WRITE(IW,9710)
            IF(CISTEP.EQ.ORMAS) THEN
               WRITE(IW,9722)
            ELSE
               WRITE(IW,9721)
            END IF
         END IF
      END IF
C
      IF(CISTEP.EQ.ALDET) CALL DETPRT(IW,NFT12,SOME)
C
      IF(CISTEP.EQ.GENCI) CALL GCIPRT(IW,NFT12,SOME)
C
C          ORMAS PRINTING REQUIRES INFO FROM DETWFN AND FCCWFN COMMONS
C
      IF(CISTEP.EQ.ORMAS) THEN
         NTOT = NCORSV + NACTDET
         NTCO = NCORSV
C
C          MODIFY MSTA TO GET RID OF CORE, WE UNDO THIS AT THE VERY END.
C
         DO II=1,NSPACE+1
            MSTA(II) = MSTA(II) - NCORSV
         ENDDO
C
         CALL VALFM(LOADFMO)
         LBOX1    = LOADFMO + 1
         LBOX2    = LBOX1   + NSPACE
         LBOX3    = LBOX2   + NSPACE
         LIBO     = LBOX3   + NSPACE
         LBST(1)  = LIBO    + NTOT
         DO II=2,NSPACE
            IDIM(II-1) = MAX(IAMA(II-1),IBMA(II-1))
            LBST(II)   = LBST(II-1) + ((MNUM(II-1)+1)*(IDIM(II-1)+1))
         ENDDO
         IDIM(NSPACE) = MAX(IAMA(NSPACE),IBMA(NSPACE))
         LAST     = LBST(NSPACE) + ((MNUM(NSPACE)+1)*(IDIM(NSPACE)+1))
         NEEDB  = LAST - LBST(1)
         NEEDO1 = LAST - LOADFM - 1
         CALL GETFM(NEEDO1)
C
         DO II=1,NSPACE
            CALL BINOM8(X(LBST(II)),MNUM(II),IDIM(II))
         ENDDO
C
C    READ ORBITALS LABELS AND GET RID OF THE CORE PART.
C
         CALL DAREAD(IDAF,IODA,X(LIBO),NTOT,262,1)
         CALL CORTRA(X(LIBO),NTOT,NTCO)
C
C    WORK OUT HOW MANY ALPHA AND BETA GROUPS THERE ARE.
C
         CALL TOTALCO(X(LBOX1),NSPACE,NADET,IAMA,IAMI,X(LBOX2),ITGA)
         CALL TOTALCO(X(LBOX1),NSPACE,NBDET,IBMA,IBMI,X(LBOX2),ITGB)
C
C    DETERMINE THE TOTAL NUMBER OF ALPHA AND BETA STRINGS.
C
         CALL RESETCO(X(LBOX1),NSPACE,NADET,IAMA,IAMI,X(LBOX2))
         IAST = 0
         DO II=1,ITGA
            CALL TOTDETG(X(LBOX1),NSPACE,LBST,X(LBST(1)),NEEDB,
     *                   MNUM,IDIM,ITOT)
            IAST = IAST + ITOT
            CALL PUSHCO(X(LBOX1),NSPACE,NADET,IAMA,IAMI,X(LBOX2),IEND)
         ENDDO
C
         CALL RESETCO(X(LBOX1),NSPACE,NBDET,IBMA,IBMI,X(LBOX2))
         IBST = 0
         DO II=1,ITGB
            CALL TOTDETG(X(LBOX1),NSPACE,LBST,X(LBST(1)),NEEDB,
     *                   MNUM,IDIM,ITOT)
            IBST = IBST + ITOT
            CALL PUSHCO(X(LBOX1),NSPACE,NBDET,IBMA,IBMI,X(LBOX2),IEND)
         ENDDO
C
C           KAB IS THE BIG STORAGE, FOR CI VECTORS
         NSYM = 2**IGPDET
         CALL VALFM(LOADFM)
         LWRK   = LOADFM + 1
         LCOA   = LWRK   + 43
         LCOB   = LCOA   + NSYM*ITGA
         LGMUL  = LCOB   + NSYM*ITGB
         LKTAB  = LGMUL  + NSYM*NSYM
         LCON   = LKTAB  + NSYM
         LCON1  = LCON   + NADET
         LCON2  = LCON1  + NADET
         LANDET = LCON2  + NADET
         LBNDET = LANDET + ITGA*NSPACE
         NAST   = LBNDET + ITGB*NSPACE
         NBST   = NAST   + ITGA+1
         LSYMA  = NBST   + ITGB+1
         LSYMB  = LSYMA  + IAST
         LGCOM  = LSYMB  + IBST
         LSPA   = LGCOM  + ITGA*ITGB
         LSPB   = LSPA   + IAST
         LDISB  = LSPB   + IBST
         LSAS   = LDISB  + NSYM*ITGB*ITGA
         LSBS   = LSAS   + (NSYM+1)*ITGA
         LSAC   = LSBS   + (NSYM+1)*ITGB
         LSBC   = LSAC   + IAST
         KAB    = LSBC   + IBST
         LAST   = KAB    + MAXP*NCIDET
         NEEDO2 = LAST - LOADFM - 1
         CALL GETFM(NEEDO2)
C
         CALL MAKTABS(IW,.FALSE.,X(LBOX1),X(LBOX2),X(LBOX3),NSPACE,
     *                NADET,NBDET,LBST,X(LBST(1)),NEEDB,
     *                IAMA,IAMI,IBMA,IBMI,MNUM,IDIM,MSTA,
     *                X(LIBO),IGPDET,KSTSYM,NSYM,
     *                NACTDET,X(LWRK),X(LKTAB),X(LGMUL),
     *                X(LCON),X(LCOA),X(LCOB),X(LANDET),X(LBNDET),
     *                X(NAST),X(NBST),X(LSYMA),X(LSYMB),X(LGCOM),
     *                MINI,MAXI,X(LSPA),X(LSPB),X(LDISB),
     *                X(LSAS),X(LSBS),X(LSAC),X(LSBC),
     *                ITGA,ITGB,IAST,IBST,NCIDET,NA1EX,NB1EX,FDIRCT)
C
         SVDSKW = DSKWRK
         DSKWRK=.FALSE.
         CALL SEQREW(NFT12)
         IF(MASWRK) READ (NFT12)
         LCIVEC = KAB
         DO IST=1,KDET
            CALL SQREAD(NFT12,X(LCIVEC),NCIDET)
            LCIVEC = LCIVEC + NCIDET
         ENDDO
         CALL SEQREW(NFT12)
         DSKWRK=SVDSKW
C
         CALL MASPRT(IW,SOME,X(KAB),X(NAST),ITGA,ITGB,X(LSYMA),
     *               IAST,IBST,X(LGCOM),X(LSBS),NSYM,X(LSBC),
     *               X(LBOX1),X(LBOX2),X(LBOX3),X(LCON1),
     *               X(LCON2),X(LKTAB))
C
         CALL RETFM(NEEDO2)
         CALL RETFM(NEEDO1)
C
C            RESTORE CORE ORBITAL COUNT BEFORE WE STOP ORMAS PRINTING
C
         DO II=1,NSPACE+1
            MSTA(II) = MSTA(II) + NCORSV
         ENDDO
      END IF
C
      IF(CISTEP.EQ.GMCCI) CALL OS_GMCCI_PRINT
C
      IF(CISTEP.EQ.GUGA  .AND.  SOME) THEN
         CALL SEQREW(NFT12)
         READ(NFT12) NSTATE,NWKS
C
         CALL SEQREW(NFT11)
         CALL RDDRT(NORBMX,NORBS,NSYM,NROWS,NWKS,LEVFRM,NEMEMX,
     *              NREFS,IEXCT,NFOCI,INTACT,NCORBS)
         IF (DROPC) NORBS = NORBS + NCORBS
C
         LEVEC  = LOADFM + 1
         LLEVIR = LEVEC  + NSTATE*NWKS
         LICASE = LLEVIR + NORBMX + 1
         LIECON = LICASE + NORBMX
         LIARC  = LIECON + NORBMX
         LAST   = LIARC  + 4*NROWS
         NEED = LAST - LOADFM
         IF(NEED.LE.NGOT) THEN
            CALL GETFM(NEED)
C
            I0=0
            DO 250 I=1,NSTATE
               CALL SQREAD(NFT12,X(LEVEC+I0),NWKS)
               I0 = I0 + NWKS
  250       CONTINUE
            CALL SEQREW(NFT12)
C
            PRTTOL = ZERO
            IF(NWKS.GT.20) PRTTOL = TOL
            IF(IPCM.EQ.1) THEN
               DO II=1,NSTATE
                  ESTATE(II)=ESTATE(II)+PC/2.0D+00
               ENDDO
            END IF
            IF(NFMOPCM.NE.0) THEN
               DO II=1,NSTATE
                  ESTATE(II)=ESTATE(II)+0.5D+00*(PC+PX)
               ENDDO
            END IF
            CALL PRCIVC(NFT11,NSTATE,NWKS,ESTATE,X(LEVEC),PRTTOL,
     *                  X(LLEVIR),X(LICASE),X(LIECON),X(LIARC),
     *                  NORBMX,4*NROWS,NROWS)
            CALL RETFM(NEED)
         ELSE
            WRITE(IW,*) 'NOT ENOUGH MEMORY TO PRINT GUGA CI EXPANSION'
         END IF
      END IF
C
C     ----- FLIP NO-S AND MO-S BACK TO THEIR PROPER SPOTS -----
C
      IF(MCNOCI) THEN
         I10=LOADFM+1
         I20=I10+N3
         I30=I20+N3
         CALL DAREAD(IDAF,IODA,X(I10),N3,15,0)
         CALL DAREAD(IDAF,IODA,X(I20),N3,19,0)
         CALL DAWRIT(IDAF,IODA,X(I10),N3,19,0)
         CALL DAWRIT(IDAF,IODA,X(I20),N3,15,0)
         RETURN
      END IF
C
C     ----- OPTIONAL EXTENDED KOOPMANS' THEOREM EIGENDATA -----
C
      IF(EKT) CALL EKTM(NORBS,ACURCY)
C
C     ----- TRANSPOSE LAGRANGIAN MATRIX -----
C     AND CHANGE ITS SHAPE TO RECTANGULAR
C
      IF (.NOT. JACOBI) THEN
         I10=LOADFM+1
         I20=I10+NLGR
         I30=I20+N3
         NEED = I30-I10
         CALL GETFM(NEED)
         CALL DAREAD(IDAF,IODA,X(I10),N3,403,0)
         DO I=1,NORBS
            DO J=1,NORB
               X(J+(I-1)*NUM+I20-1)=X(I+(J-1)*NORBS+I10-1)
            ENDDO
         ENDDO
         CALL DAWRIT(IDAF,IODA,X(I20),N3,403,0)
         CALL RETFM(NEED)
      ENDIF
C
C     ----- GENERATE NATURAL ORBITALS AND DENSITY IN AO BASIS -----
C
      LVEC  = LOADFM+1
      LQ    = LVEC  + N3
      LGIJ  = LQ    + N3
      LGMN  = LGIJ  + N2
      LOCC  = LGMN  + N2
      LWRK  = LOCC  + N1
      LIWRK = LWRK  + N1
      LSCR  = LIWRK + N1
      LSQ   = LSCR  + N1*8
      LAST  = LSQ   + N3
      NEED  = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
C                TRANSFER DM1 FROM NFT15 TO DAF RECORD 68
C                MIMICKS NTNONE SKIPPED BY PARALLEL DIRECT HESS
C
      IF (GOPARR.AND.FULLNR) CALL PNTNWT(NFT15,X(LGIJ),L2)
C
      CALL NTNDEN(X(LVEC),X(LQ),X(LGIJ),X(LGMN),X(LOCC),X(LWRK),
     *            X(LIWRK),X(LSCR),X(LSQ),N0,N1,N2,N3,L1,L2,
     *            NCORBS,SOME)
      CALL RETFM(NEED)
C
C     ----- MAKE CANONICAL ORBITALS -----
C     OR, IN ANY CASE, PRINT OUT THE NO-S AND MO-S
C
C     CANONICALIZATION TAKES SOME TIME, SO WE CAN SKIP IT DURING
C     THE COURSE OF A GEOMETRY SEARCH OR NUMERICAL HESSIAN.
C     CANONICALIZATION CAN DAMAGE ORBITAL SYMMETRY OF CORES,
C     ESPECIALLY IN NON-ABELIAN GROUPS.
C     NOTE THAT MR-PT OR MR-CI METHODS ARE BASED ON CANONICALIZATION
C     HAVING SEPARATED CORE FROM VALENCE, SO DO NOT SKIP THIS,
C     UNLESS THE USER'S INPUT HAS REQUESTED THAT.
C
      CANSAV = CANONC
      IF(RUNTYP.EQ.OPT  .OR.  RUNTYP.EQ.SAD  .OR.  IHESSM.EQ.1)
     *          CANONC = .FALSE.
      IF(MPLEVL.GT.0)    CANONC=CANSAV
      IF(CITYP.NE.RNONE) CANONC=CANSAV
C
      NVAL  = NORBS-NCORBS
      NVIRT = NORB-NORBS
      CALL CANMCV(NPFLG8,CANONC,NPUNCH,NORB,NCORBS,NVAL,NVIRT,NORBS)
      CANONC = CANSAV
C
C     ----- SWITCH -MCHF- AND -NO- ORBS IF -NOCI- = .TRUE. -----
C     THIS PERMITS AN ENSUING CI TO BE BASED ON NATURAL ORBITALS.
C
      IF(NOCI) THEN
         I10=LOADFM+1
         I20=I10+NUM*NUM
         I30=I20+NUM*NUM
         CALL DAREAD(IDAF,IODA,X(I10),N3,15,0)
         CALL DAREAD(IDAF,IODA,X(I20),N3,19,0)
         CALL DAWRIT(IDAF,IODA,X(I10),N3,19,0)
         CALL DAWRIT(IDAF,IODA,X(I20),N3,15,0)
      END IF
      RETURN
C
 9000 FORMAT(/5X,32(1H-),5X,29(1H-)/
     *        5X,'MCHF NEWTON-RAPHSON OPTIMIZATION',5X,
     *           'CODE WRITTEN BY MICHEL DUPUIS'/
     *        5X,32(1H-),5X,29(1H-))
 9010 FORMAT(1X,'MEMORY DEMAND FOR 2 ELECTRON NTN-RPH CONTRIBUTIONS'/
     *       1X,'     METHOD=DM2     NEEDS',I10,' WORDS'/
     *       1X,'     METHOD=TEI     NEEDS',I10,' WORDS'/
     *       1X,'TOTAL AVAILABLE MEMORY IS',I10,' WORDS')
C
 9100 FORMAT(1X,'CHOOSING DM2 DRIVEN METHOD...')
 9110 FORMAT(1X,'NOT ENOUGH CORE FOR METHOD=DM2 NTN-RPH MATRICES.'/
     *       1X,'NEED,NGOT = ',2I10,' WORDS.')
 9120 FORMAT(1X,'FORMING -LGR- + -COU- + -EXC- FROM',
     *          ' (IJ//KL) INTEGRALS')
 9130 FORMAT(1X,'FORMING -LGR- + -COU- + -EXC- FROM',
     *          ' (AJ//KL) INTEGRALS')
 9140 FORMAT(1X,'FORMING         -COU-         FROM',
     *          ' (AB//KL) INTEGRALS')
 9150 FORMAT(1X,'FORMING                 -EXC- FROM',
     *          ' (AJ//BL) INTEGRALS')
C
 9200 FORMAT(1X,'CHOOSING INTEGRAL DRIVEN METHOD...')
 9210 FORMAT(1X,'NOT ENOUGH CORE FOR METHOD=TEI NTN-RPH MATRICES.'/
     *       1X,'NEED,NGOT = ',2I10,' WORDS.')
C
 9500 FORMAT(1X,'..... DONE WITH 2 ELECTRON CONTRIBUTIONS .....')
 9510 FORMAT(1X,'NOT ENOUGH CORE TO BUILD 1E-LAGRANGIAN MATRIX.'/
     *       1X,'NGOT,NEED = ',2I10,' WORDS.')
 9520 FORMAT(1X,'..... DONE WITH 1 ELECTRON CONTRIBUTIONS .....')
 9530 FORMAT(1X,'ITER = ',I4,' MICRO = ',I4,
     *          ' ENERGY = ',F19.9,' DEMAX = ',F16.9)
 9540 FORMAT(1X,'-MROT- IS FOUND GREATER THAN -NROT-.',2I5)
 9550 FORMAT(1X,'NOT ENOUGH CORE TO SET UP NEWTON-RAPHSON MATRICES.'/
     *       1X,'NGOT,LAST = ',2I10,' WORDS.')
 9560 FORMAT(1X,'..... DONE SOLVING NEWTON-RAPHSON EQUATIONS .....')
 9570 FORMAT(1X,'NOT ENOUGH CORE TO BUILD ORBITAL ROTATION MATRIX.'/
     *       1X,'NGOT,NEED = ',2I10,' WORDS.')
 9580 FORMAT(1X,'..... DONE WITH NTN-RPH ORBITAL ROTATIONS ..... ')
 9590 FORMAT(1X,'MAXIMUM MEMORY USED BY NTN-RPH CODE WAS',I10,' WORDS.')
C
 9700 FORMAT(/1X,'-MCCI- BASED ON NATURAL ORBITALS'/1X,32(1H-))
 9710 FORMAT(/1X,'-MCCI- BASED ON OPTIMIZED ORBITALS'/1X,34(1H-))
 9721 FORMAT(/1X,'PLEASE NOTE: IF THE ACTIVE ORBITALS ARE',
     *           ' CANONICALIZED BELOW,'/
     *        1X,'THE FOLLOWING CI EXPANSION COEFFICIENTS AND THE',
     *           ' DENSITY DO NOT'/
     *        1X,'CORRESPOND TO THE PRINTED ORBITALS.  THE PRINTED',
     *           ' EXPANSIONS MATCH'/
     *        1X,'THE ORBITALS USED DURING THE LAST ITERATION.  IF',
     *           ' YOU WISH TO SEE'/
     *        1X,'CI EXPANSIONS BASED ON THE CANONICAL (OR NATURAL)',
     *           ' ORBITALS, YOU'/
     *        1X,'MUST RUN A CI CALCULATION WITH THAT ORBITAL CHOICE',
     *           ' READ IN $VEC.')
 9722 FORMAT(/1X,'PLEASE NOTE: ORMAS IS NOT INVARIANT TO ORBITAL',
     *           ' ROTATIONS BETWEEN THE'/
     *        1X,'DIFFERENT ACTIVE SPACES.  THE CI EXPANSION JUST',
     *           ' BELOW IS FOR THE'/
     *        1X,'OPTIMIZED ORBITALS USED IN THE LAST ITERATION. ',
     *           ' THE NATURAL ORBITALS'/
     *        1X,'ARE GENERATED BELOW, BUT THE NO''S ROTATE BETWEEN',
     *           ' THE ACTIVE SPACES,'/
     *        1X,'SO THE NATURAL ORBITALS DO NOT REPRODUCE THE',
     *           ' ORMAS WAVEFUNCTION.')
 9730 FORMAT(1X,'ACTIVE ORBITAL FOCK MATRIX CONSTRUCTION NEEDS ',I10,
     *          ' WORDS',/
     *       1X,'AVAILABLE MEMORY IS ',I10,' WORDS')
 9740 FORMAT(1X,'FORMING FOCK OPERATOR OVER ACTIVE ORBITAL DENSITY,',
     *          ' MEMORY=',I10)
 9900 FORMAT(5X,'PARALLEL VERSION WRITTEN BY G. FLETCHER (ELORET)'/
     *       5X,'NASA ADVANCED SUPERCOMPUTING DIVISION,',
     *          ' AMES RESEARCH CENTER'/)
      END
C*MODULE MCSCF   *DECK NTNCAN
      SUBROUTINE NTNCAN(V,VOLD,FMO,F,S,Q,U,OCC,EIG,SCR,IWRK,
     *                  NORB,NCORBS,NVAL,NVIRT,NORBS,L0,L1,L2,L3,
     *                  ETOT,ENUCR,SOME,DBG,NPUNCH)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION METHOD
C
      PARAMETER (MXATM=2000, MXAO=8192, MXNORO=250)
C
      LOGICAL FORS,SOME,DBG,GOPARR,DSKWRK,MASWRK,JACOBI,
     *        CANONC,FCORE,FORSNR,NOCI,EKT,LINSER,FDIRCT,QCORR
C
      DIMENSION V(L1,L1),VOLD(L1,L1),FMO(L2),F(L2),S(L2),Q(L3),U(L1,L1),
     *          OCC(L1),EIG(L1),SCR(L1,8),IWRK(L1)
      DIMENSION TIMSTR(3)
C
      COMMON /FCCWFN/ NSPACE,MSTA(51),MNUM(51),MINI(51),MAXI(51),
     *                IAMI(51),IAMA(51),IBMI(51),IBMA(51),IDIM(51),
     *                LBST(51),NREF0,FDIRCT,QCORR,C0SQ
      COMMON /GUGWFN/ NGFZC,NGMCC,NGDOC,NGAOS,NGBOS,NGALP,NGVAL,
     *                NGEXT,NGFZV,IFORS,IEXCIT,ICICI,NOIRR
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /JACOBI/ JACOBI,NJAOR,ELAST,ISTAT
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,
     *                NNORB,NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORSNR,NOCI,EKT,LINSER
      COMMON /MEXOPT/ MEXSKPGES,MEXSTATE
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PRPOPT/ ILOCAL
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00)
C
      DATA VECWRD,VEC1,VEC2,ENDWRD
     *     /8H $VEC   ,8H $VEC1  ,8H $VEC2  ,8H $END   /
      DATA GUGA,ALDET,GENCI,ORMAS,GMCCI
     *     /8HGUGA    ,8HALDET   ,8HGENCI   ,8HORMAS   ,8HGMCCI   /
      DATA AMEX/8HMEX     /
C
C          NCORBS,NVAL,NVIRT ARE SIZES OF CORE,ACTIVE,EXTERNAL SPACES
C          NORBS IS THE NUMBER OF OCCUPIED ORBITALS, WHILE
C          NORB IS THE TOTAL NUMBER OF MOS USED IN THE VARIATION SPACE
C            (USUALLY NORB WILL EQUAL L0, UNLESS RESET IN $MCSCF)
C
      IF((NORB   .NE.  (NCORBS+NVAL+NVIRT))  .OR.
     *    NORBS  .NE.  (NCORBS+NVAL)) THEN
         IF(MASWRK) WRITE(IW,9010) NORB,NORBS,NCORBS,NVAL,NVIRT
         CALL ABRT
      END IF
C
C         THIS ROUTINE MAY ALTER CORE AND VIRTUAL ORBITALS,
C         AND IT MAY CHANGE THE PHASE OF ACTIVE ORBITALS.
C         SINCE THE CI VECTOR WILL NOT BE UPDATED TO REFLECT
C         THESE POSSIBLE PHASE CHANGES, WE MUST SAVE AN ABSOLUTELY
C         PRISTINE SET OF THE FINAL ORBITALS IN CASE THEY'RE NEEDED.
C
      CALL DAREAD(IDAF,IODA,V,L3, 15,0)
      CALL DAWRIT(IDAF,IODA,V,L3,325,0)
C
C         IS THIS MCSCF OF THE FULL ACTIVE SPACE TYPE?
C
      IF(CISTEP.EQ.GUGA) THEN
         NACTE = 2*NGDOC + NGAOS + NGBOS + NGALP
         NHOLE =           NGAOS + NGBOS + NGALP + 2*NGVAL
         IEXCT = MIN(NACTE,NHOLE)
         FORS = IEXCT.LE.IEXCIT  .AND.  NGEXT.EQ.0
      END IF
      IF(CISTEP.EQ.ALDET) FORS=.TRUE.
      IF(CISTEP.EQ.ORMAS) THEN
         FORS=.FALSE.
         IF(NSPACE.EQ.1) FORS=.TRUE.
      END IF
      IF(CISTEP.EQ.GENCI) FORS=.FALSE.
      IF(CISTEP.EQ.GMCCI) FORS=.FALSE.
C
C         IF A GRADIENT IS TO BE COMPUTED, THE ACTIVE SPACE MUST
C         NOT BE ROTATED HERE, OR THE LAGRANGIAN WILL BE INCORRECT.
C         THE LAGRANGIAN IS DEFINED FROM THE TWO PARTICLE DENSITY,
C         AND THIS MAY BE WHY THE CODE BELOW DOESN'T ROTATE THE
C         LAGRANGIAN CORRECTLY???  SOMEDAY THIS SHOULD BE REVISITED.
C
      CALL DERCHK(NDER)
      IF(NDER.GT.0) FORS=.FALSE.
C
C         MRMP/MCQDPT HAS NO ANALYTIC GRADIENT, AND IS NECESSARILY
C         BASED ON FCI.  IT IS EFFICIENT TO CANONICALIZE THE ACTIVE
C         SPACE HERE, TO ALLOW THE FULL CI REFERENCE SPACE MRMP/MCQDPT
C         CODES TO SKIP CANONICALIZING (IT SAVES A TRANSFORMATION).
C         NOTE THAT -CANONC- INPUT MAY STILL SKIP ALL CANONICALIZATIONS.
C
      IF(CISTEP.EQ.ALDET  .AND.  MPLEVL.EQ.2) FORS=.TRUE.
C
C         LOCALIZED ORBITAL ANALYSIS USES THE MO DENSITY ON RECORD 68,
C         AND THIS ROUTINE DOES NOT UPDATE THIS DAF RECORD AT PRESENT.
C         ACTUALLY, THIS MIGHT BE RELEVANT TO THE GRADIENT PROBLEM
C         MENTIONED IMMEDIATELY ABOVE.
C
      IF(ILOCAL.NE.0) FORS=.FALSE.
C
      IF(CANONC) THEN
         IF(FORS) THEN
            IF(SOME) WRITE(IW,9020) NCORBS,NVAL,NVIRT
         ELSE
            IF(SOME) WRITE(IW,9030) NCORBS,NVIRT
         END IF
      ELSE
         IF(SOME) WRITE(IW,9040)
         CALL DAREAD(IDAF,IODA,  V,L3,19,0)
         CALL DAREAD(IDAF,IODA,OCC,L1,21,0)
         GO TO 600
      END IF
C
C     ----- TRANSFORM THE FOCK OPERATOR TO MCSCF MO BASIS -----
C
      CALL DAREAD(IDAF,IODA,V,L3,15,0)
      CALL TFTRI(FMO,F,V,SCR,NORB,L1,L1)
      IF(DBG) THEN
         IF(MASWRK) WRITE(IW,9050)
         CALL PRTRI(FMO,NORB)
      END IF
C
C     ----- ZERO BLOCKS CONNECTING CORE, ACTIVE, EXTERNAL SPACES -----
C     THIS WILL PREVENT ANY SUBSEQUENT MIXING OF THESE ORBITAL SPACES
C
      DO I=NCORBS+1,NCORBS+NVAL
         IROW = (I*I-I)/2
         DO J=1,NCORBS
            IJ = IROW+J
            FMO(IJ) = ZERO
         ENDDO
      ENDDO
      DO I=NCORBS+NVAL+1,NORB
         IROW = (I*I-I)/2
         DO J=1,NCORBS+NVAL
            IJ = IROW+J
            FMO(IJ) = ZERO
         ENDDO
      ENDDO
C
C     ----- TRANSFORM BLOCKED FOCK MATRIX BACK TO AO SPACE -----
C
      CALL DAREAD(IDAF,IODA,S,L2,12,0)
      CALL TFTRIB(F,FMO,S,V,U,SCR,L0,L1,L2,L3)
C
C     ----- DIAGONALIZE IN THE -Q- MATRIX SPACE -----
C
      CALL DAREAD(IDAF,IODA,Q,L3,45,0)
      CALL TFTRI(FMO,F,Q,SCR,L0,L1,L1)
      CALL SYMDIA(FMO,V,EIG,SCR,IWRK,L0,L2,L1)
      CALL TFSQB(V,Q,SCR,L0,L1,L1)
C
C        COPY THE ORIGINAL MCSCF ORBITALS TO UPDATE THE LAGRANGIAN LATER
C
      CALL DAREAD(IDAF,IODA,U,L3,15,0)
      CALL DCOPY(L3,U,1,VOLD,1)
C
C        MAKE SURE WE PRESERVE ORIGINAL CORE/ACTIVE/EXTERNAL ORDERING.
C        THEN WITHIN EACH ORBITAL SPACE, WE SORT BY EIGENVALUE.
C
      CALL OVLSEL(V,EIG,VOLD,S,U,SCR,IWRK,SCR,L0,L1,L2)
C
      DO K=1,NCORBS-1
         KMIN=K
         RMIN=EIG(K)
         DO I=K+1,NCORBS
            IF(EIG(I).LT.RMIN) THEN
               KMIN=I
               RMIN=EIG(I)
            END IF
         END DO
         IF(KMIN.NE.K) THEN
            CALL DSWAP( 1,EIG(K),1,EIG(KMIN),1)
            CALL DSWAP(L1,V(1,K),1,V(1,KMIN),1)
         END IF
      ENDDO
C
      IF(FORS) THEN
      DO K=NCORBS+1,NORBS-1
         KMIN=K
         RMIN=EIG(K)
         DO I=K+1,NORBS
            IF(EIG(I).LT.RMIN) THEN
               KMIN=I
               RMIN=EIG(I)
            END IF
         END DO
         IF(KMIN.NE.K) THEN
            CALL DSWAP( 1,EIG(K),1,EIG(KMIN),1)
            CALL DSWAP(L1,V(1,K),1,V(1,KMIN),1)
         END IF
      END DO
      END IF
C
      DO K=NORBS+1,NORBS+NVIRT-1
         KMIN=K
         RMIN=EIG(K)
         DO I=K+1,NORBS+NVIRT
            IF(EIG(I).LT.RMIN) THEN
               KMIN=I
               RMIN=EIG(I)
            END IF
         END DO
         IF(KMIN.NE.K) THEN
            CALL DSWAP( 1,EIG(K),1,EIG(KMIN),1)
            CALL DSWAP(L1,V(1,K),1,V(1,KMIN),1)
         END IF
      END DO
C
C        NOW APPLY A CONSISTENT CHOICE FOR THE PHASE FACTOR.
C        NOTE THAT IF WE REPLACE ANY ORBITALS BY THE ORIGINALS,
C        THE REPLACEMENTS MUST NOT HAVE ANY PHASES CHANGED!
C
      CALL STFASE(V,L1,L1,L0)
C
C        FOR FROZEN CORE MCSCF, REPLACE CANONICAL CORE BY THE ORIGINALS
C
      IF(FCORE) THEN
         DO K=1,NCORBS
            CALL DCOPY(L1,VOLD(1,K),1,V(1,K),1)
            EIG(K) = ZERO
         ENDDO
      END IF
C
C        IN CASE THIS IS NOT A FULL ACTIVE SPACE, RE-INSERT THE
C        OPTIMAL MOS, RATHER THAN USING CANONICAL ACTIVE ORBITALS.
C        PRIOR TO NOVEMBER 2001, GAMESS ALWAYS RE-INSERTED HERE.
C
      IF(.NOT.FORS) THEN
         DO K=NCORBS+1,NCORBS+NVAL
            CALL DCOPY(L1,VOLD(1,K),1,V(1,K),1)
            EIG(K) = ZERO
         ENDDO
      END IF
C
      CALL DAWRIT(IDAF,IODA,  V,L3,15,0)
      CALL DAWRIT(IDAF,IODA,EIG,L1,17,0)
      IF(DBG) THEN
         IF(MASWRK) WRITE(IW,9060)
         CALL PREV(V,EIG,L0,L0,L1)
      END IF
C
C        INSERT THE ORIGINAL ACTIVE NATURAL ORBITALS AND OCCUPANCIES
C
      CALL DAREAD(IDAF,IODA,  U,L3,19,0)
      CALL DAREAD(IDAF,IODA,OCC,L1,21,0)
      DO K=NCORBS+1,NCORBS+NVAL
         EIG(K) = OCC(K)
         DO I=1,L1
            V(I,K) = U(I,K)
         ENDDO
      ENDDO
      CALL DCOPY(L1,EIG,1,OCC,1)
      CALL DAWRIT(IDAF,IODA,  V,L3,19,0)
C----      CALL DAWRIT(IDAF,IODA,OCC,L1,21,0)
C
C     ----- PRINT MCSCF NATURAL ORBITAL SET -----
C
  600 CONTINUE
      IF(SOME) THEN
         WRITE(IW,9070)
         CALL PREV(V,OCC,NORBS,L1,L1)
      END IF
      IF(MASWRK  .AND.  NPUNCH.GT.0) THEN
         CALL TMDATE(TIMSTR)
         WRITE(IP,8000) TIMSTR
         WRITE(IP,8020) TITLE,ETOT,ENUCR
         IF(RUNTYP.NE.AMEX) THEN
            WRITE(IP,8030) VECWRD
         ELSE
            IF(MEXSTATE.EQ.1) WRITE(IP,8030) VEC1
            IF(MEXSTATE.EQ.2) WRITE(IP,8030) VEC2
         END IF
         CALL PUSQL(V,NORBS,L1,L1)
         WRITE(IP,8030) ENDWRD
      END IF
C
C     ----- ASSIGN SYMMETRY LABELS TO -MCHF- ORBITALS -----
C
      CALL DAREAD(IDAF,IODA,  S,L2,12,0)
      CALL DAREAD(IDAF,IODA,  V,L3,15,0)
      CALL DAREAD(IDAF,IODA,  Q,L3,45,0)
      CALL SYMMOS(IWRK,Q,S,V,SCR,L0,L1,L0,L1)
      CALL DAWRIT(IDAF,IODA,IWRK,L1,255,1)
C
C     ----- PRINT -MCHF- ORBITALS -----
C
      IF(CANONC) THEN
         CALL DAREAD(IDAF,IODA,EIG,L1,17,0)
      ELSE
         CALL VCLR(EIG,1,L1)
      END IF
      IF(SOME) THEN
         WRITE(IW,9080)
         CALL PREVS(V,EIG,IWRK,L0,L1,L1,1)
      END IF
C
      IF(MASWRK  .AND.  NPUNCH.GT.0) THEN
         WRITE(IP,8010) TIMSTR
         WRITE(IP,8020) TITLE,ETOT,ENUCR
         IF(RUNTYP.NE.AMEX) THEN
            WRITE(IP,8030) VECWRD
         ELSE
            IF(MEXSTATE.EQ.1) WRITE(IP,8030) VEC1
            IF(MEXSTATE.EQ.2) WRITE(IP,8030) VEC2
         END IF
         NPUN = L0
         IF(NPUNCH.EQ.1) NPUN=NORBS
         CALL PUSQL(V,NPUN,L1,L1)
         WRITE(IP,8030) ENDWRD
      END IF
      IF(.NOT.CANONC) RETURN
      IF(JACOBI) RETURN
C
C     -- UPDATE THE LAGRANGIAN MATRIX TO REFLECT THE ORBITAL CHANGES --
C     SINCE WE DON'T HAVE A SIMPLE ROTATION MATRIX FROM THE OLD TO NEW
C     ORBITALS, TRANSFORM BACK TO THE AO BASIS, AND THEN FORWARD INTO
C     THE NEW MO BASIS.  RELOADING OF THE ORIGINAL CANONICAL ORBITALS
C     MAKES SURE THAT ANY SIGN CHANGES IN -SYMMOS- ARE NOT CARRIED ON.
C
C     THERE IS GARBAGE IN THE ANY-ORBITAL/VIRTUAL BLOCK AT THIS POINT,
C     BUT THESE ELEMENTS ARE OF COURSE IDENTICALLY ZERO BY DEFINITION.
C
      CALL DAREAD(IDAF,IODA,V,L3,15,0)
      CALL DAREAD(IDAF,IODA,U,L3,403,0)
      CALL VCLR(U(1,NORBS+1),1,L1*(NORB-NORBS))
C
C        TRANSFORM LAGRANGIAN TO AO BASIS, E' = VOLD * E * VOLD-DAGGER
C
      CALL MRARBR(VOLD,L1,L1,NORB,U,L1,NORB,Q,L1)
      CALL MRARTR(Q,L1,L1,NORB,VOLD,L1,L1,U,L1)
C
C        SV PRODUCT WILL BE STORED IN Q
C
      CALL DAREAD(IDAF,IODA,S,L2,12,0)
      CALL MTARBR(S,L1,V,NORB,Q,L1,1)
C
C        TRANSFORM TO CANONICAL MO SPACE, E = (SV)-DAGGER * E' * (SV)
C        NOTE THAT THE ORIGINAL ORBITALS ARE USED AS A WORK SPACE.
C
      CALL MRTRBR(Q,L1,L1,NORB,U,L1,L1,VOLD,L1)
      CALL MRARBR(VOLD,L1,NORB,L1,Q,L1,NORB,U,L1)
      CALL DAWRIT(IDAF,IODA,U,L3,403,0)
      RETURN
C
 8000 FORMAT('--- NATURAL ORBITALS OF MCSCF --- GENERATED AT ',3A8)
 8010 FORMAT('--- OPTIMIZED MCSCF MO-S --- GENERATED AT ',3A8)
 8020 FORMAT(10A8/'E(MCSCF)=',F20.10,', E(NUC)=',F20.10)
 8030 FORMAT(A8)
C
 9010 FORMAT(1X,'*** ERROR IN NTNCAN, CHECK CALLING PARAMETERS:'/
     *       1X,'NUM. OCCUPIED ORBS -NORBS- IS',I5,
     *          ' TOTAL NUM. ORBS -NORB- IS',I5/
     *       1X,'BUT NUM. CORE, ACTIVE, VIRTUAL=',3I5)
 9020 FORMAT(1X,I5,' FILLED,',I5,' ACTIVE, AND',I5,' VIRTUAL ORBITALS',
     *           '  WILL BE CANONICALIZED')
 9030 FORMAT(1X,'ONLY THE',I5,' CORE AND',I5,' VIRTUAL ORBITALS',
     *          ' WILL BE CANONICALIZED')
 9040 FORMAT(/1X,'NO ORBITALS WILL BE CANONICALIZED')
 9050 FORMAT(/1X,'THE STANDARD FOCK OPERATOR IN THE MO BASIS IS')
 9060 FORMAT(/1X,'THE CORE/ACTIVE/VIRTUAL CANONICAL ORBITALS ARE')
 9070 FORMAT(/10X,22(1H-)/10X,'MCSCF NATURAL ORBITALS'/
     *        10X,22(1H-))
 9080 FORMAT(/10X,24(1H-)/10X,'MCSCF OPTIMIZED ORBITALS'/
     *        10X,24(1H-))
      END
C*MODULE MCSCF   *DECK NTNDEN
      SUBROUTINE NTNDEN(VEC,Q,GIJ,GMN,OCC,WRK,IWRK,SCR,SQ,
     *                  N0,N1,N2,N3,L1,L2,NCORBS,SOME)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SOME,GOPARR,DSKWRK,MASWRK
C
C           HERE -L- MEANS OCCUPIED AND -N- ENTIRE ORBITAL SPACE.
C           STORAGE -GIJ- IS USED AS BOTH -L2- AND -N2- SIZES
      DIMENSION VEC(N1,N1),Q(N1,N1),GIJ(*),GMN(N2),OCC(N1),
     *          WRK(N1),IWRK(N1),SCR(N1,8),SQ(N3)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,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
      DATA RNONE/8HNONE    /
C
C     ----- GET NATURAL ORBITALS OF MCSCF WAVEFUNCTION -----
C     CAUTION: -L- INTEGERS ARE OCCUPIED MO SPACE, AND -N- THE AO SPACE.
C
C        BEGIN BY TRANSFORMING THE DENSITY MATRIX TO THE AO BASIS.
C        THIS IS SAVED FOR SUBSEQUENT PROPERTY EVALUATION.
C
      CALL DAREAD(IDAF,IODA,VEC,N3,15,0)
      CALL DAREAD(IDAF,IODA,GIJ,L2,68,0)
      CALL TRPOSQ(VEC,N1)
      CALL TFTRI(GMN,GIJ,VEC,WRK,N1,L1,N1)
      CALL DAWRIT(IDAF,IODA,GMN,N2,16,0)
C
C        TRANSFORM THE AO BASIS DENSITY MATRIX TO Q-VECTOR SPACE,
C        USING THE TRANSFORMATION MATRIX S*Q, AND DIAGONALIZE IT.
C        THE SIGN CHANGE IS TO SORT EIGENVALUES INTO ASCENDING ORDER.
C
      CALL DAREAD(IDAF,IODA,GIJ,N2,12,0)
      CALL DAREAD(IDAF,IODA,Q  ,N3,45,0)
      CALL MTARBR(GIJ,N1,Q,N0,SQ,N1,1)
      CALL TFTRI(GIJ,GMN,SQ,WRK,N0,N1,N1)
      DO 100 I=1,(N0*N0+N0)/2
         GIJ(I)=-GIJ(I)
  100 CONTINUE
      KDSAVE = KDIAG
      KDIAG = 3
      CALL SYMDIA(GIJ,VEC,OCC,SCR,IWRK,N0,N2,N1)
      KDIAG = KDSAVE
      CALL TFSQB(VEC,Q,WRK,N0,N1,N1)
      DO 120 I=1,N0
         OCC(I)=-OCC(I)
  120 CONTINUE
C
C        THE NATURAL ORBITALS ARE PRINTED LATER,
C        AFTER CANONICALIZATION IN THE CORE SPACE.
C
      CALL DAWRIT(IDAF,IODA,OCC,N1,21,0)
      CALL DAWRIT(IDAF,IODA,VEC,N3,19,0)
C
      NVAL = L1-NCORBS
      MINI = NCORBS+1
      MAXI = NCORBS+NVAL
C
C     ----- PRINT DENSITY OVER ACTIVE MO'S -----
C
      CALL DAREAD(IDAF,IODA,GIJ,L2,68,0)
      IF(SOME) THEN
         IJ = (MINI*MINI+MINI)/2
         LOC = 1
         DO 650 I=MINI,MAXI
            NCOL = I - NCORBS
            CALL DCOPY(NCOL,GIJ(IJ),1,GMN(LOC),1)
            IJ = IJ + I
            LOC = LOC + NCOL
  650    CONTINUE
         WRITE(IW,9020)
         CALL PRTRI(GMN,NVAL)
      END IF
C
C        CHECK TO SEE IF THE ACTIVE SPACE HAS 2 OR 0 OCCUPANCIES,
C        AND PRINT A WARNING IF IT DOES.
C
      NZERO=0
      NDOCC=0
      DO I=MINI,MAXI
         IF(ABS(OCC(I))        .LT.1.0D-04) NZERO=NZERO+1
         IF(ABS(OCC(I)-2.0D+00).LT.1.0D-04) NDOCC=NDOCC+1
      ENDDO
      IF((NDOCC.GT.0  .OR.  NZERO.GT.0)  .AND.  MASWRK) THEN
         WRITE(IW,9030)
         IF(NDOCC.GT.0) WRITE(IW,9040) NDOCC
         IF(NZERO.GT.0) WRITE(IW,9050) NZERO
         WRITE(IW,9060)
         IF(MPLEVL.GT.0  .OR.  CITYP.NE.RNONE) THEN
            WRITE(IW,9030)
            WRITE(IW,9070)
         ENDIF
      END IF
      RETURN
C
 9020 FORMAT(/1X,'DENSITY MATRIX OVER ACTIVE MO-S')
 9030 FORMAT(/1X,'* * * WARNING * * *')
 9040 FORMAT(1X,'YOUR NATURAL ORBITALS HAVE',I4,
     *          ' DOUBLY OCCUPIED ACTIVE ORBITALS.')
 9050 FORMAT(1X,'YOUR NATURAL ORBITALS HAVE',I4,
     *          ' UNOCCUPIED ACTIVE ORBITALS.')
 9060 FORMAT(1X,'AT THE VERY LEAST, THIS IS A WASTE OF COMPUTER TIME,'/
     * 1X,'BUT IT MAY REVEAL A FLAW IN YOUR MCSCF COMPUTATION:'/
     * 3X,'CHECK YOUR FORMULATION OF THE ACTIVE SPACE,'/
     * 3X,'AS WELL AS MECHANICAL THINGS LIKE STARTING ORBITALS.')
 9070 FORMAT(1X,'THIS RUN REQUESTS EITHER MPLEVL=2 OR CITYP=XXX ENERGY',
     *          ' CORRECTIONS.'/
     * 1X,'THE INCLUSION OF FILLED AND/OR EMPTY ORBITALS IN THE ACTIVE'/
     * 1X,'SPACE WILL RESULT IN INCORRECT DEFAULT ORBITAL',
     *    ' CANONICALIZATION'/
     * 1X,'OF THE CORE AND/OR VIRTUAL SPACE.'/
     * 1X,'THEREFORE, THE DYNAMIC CORRELATION ENERGY COMPUTED FOR',
     *    ' ANY FILLED'/
     * 1X,'MCSCF VALENCE ORBITALS WILL BE INCORRECT, BECAUSE THESE',
     *    ' WILL NOT'/
     * 1X,'HAVE BEEN CORRECTLY SEPARATED FROM THE CHEMICAL CORE.')
      END
C*MODULE MCSCF   *DECK NTNDMP
      SUBROUTINE NTNDMP(DE,DEP,DEAVG,DAMP)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00, FOUR=4.0D+00,
     *           PT25=2.5D-01, PT5=5.0D-01, PT2=0.2D-01,
     *           FAC=1.6D+01)
C
      IF( DE.GT.ZERO) GO TO 400
      IF(DEP.GT.ZERO) GO TO 300
      IF( DE.GT. DEP) GO TO 200
C
C     ----- DE < 0. , DEP < 0. , DE < DEP -----
C
      IF( ABS(DE).LT.TWO*DEAVG) GO TO 110
      DAMP=FAC* MAX(DAMP,DEAVG)
      RETURN
  110 IF( ABS(DE).GT.PT5*DEAVG) RETURN
      DAMP=DAMP/FAC
      RETURN
  200 CONTINUE
C
C     ----- DE < 0. , DEP < 0. , DE > DEP -----
C
      IF(DE.GT.PT25*DEP) GO TO 210
      DAMP=(DE/DEP)**2* MAX(DAMP,DEAVG)
      RETURN
  210 DAMP=DAMP/FAC
      RETURN
  300 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) RETURN
      DAMP=DAMP/FAC
      RETURN
  400 CONTINUE
      IF(DEP.GT.ZERO) GO TO 500
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) RETURN
      DAMP=DAMP/FAC
      RETURN
  500 CONTINUE
C
C     ----- DE > 0. , DEP > 0. -----
C
      DAMP=FOUR* MAX(DAMP,DEAVG)
      IF(DE.LT.FOUR*DEP) GO TO 510
      DAMP=FAC* MAX(DAMP,DEAVG)
      RETURN
  510 IF(DE.GT.PT25*DEP) GO TO 520
      DAMP=DAMP/FAC
      RETURN
  520 DAMP=(DE/DEP)**2* MAX(DAMP,DEAVG)
      RETURN
      END
C*MODULE MCSCF   *DECK NTNFCK
      SUBROUTINE NTNFCK(D,F,XX,IX,NINTMX,A,IA,NUM,NOSQUR)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL IANDJ,KANDL,DIRTRF
C
      DIMENSION D(*),F(*),IA(NUM),A(*)
      DIMENSION XX(NINTMX),IX(NINTMX)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /PCKLAB/ LABSIZ
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
      COMMON /INT2IC/ NINTIC,ININTIC,NXXIC,LBUFPIC,LIXIC,LABSIX,NINTIX
C
      PARAMETER (ZERO=0.0D+00, PT5=0.5D+00)
C
C     ----- FORM FOCK OPERATOR -----
C     THE CALCULATION IS DONE IN THE AO BASIS, -D- IS THE TOTAL
C     ONE PARTICLE DENSITY, IJKO IS THE ORDERED AO INTEGRAL FILE.
C
      N2=(NUM*NUM+NUM)/2
      CALL VCLR(F,1,N2)
      CALL DAREAD(IDAF,IODA,D,N2,16,0)
C
C     ----- INTEGRALS IN TRIANGULAR ORDER -----
C
      IF(NOSQUR.EQ.1) GO TO 1000
      CALL SEQREW(IJKO)
C     FIRST, PROCESS IN CORE INTEGRALS
      NPACKIC=0
  100 CONTINUE
      IF(NPACKIC.EQ.0.AND.NINTIC.NE.0) THEN
        NXX=NXXIC
        NINT=NXX
      ELSE
C       IF ALL IN CORE BAIL OUT
        IF(NINTIC.NE.0.AND.NXXIC.LT.NINTIC) GO TO 400
        CALL PREAD(IJKO,XX(NINTIC+1),IX(ININTIC+1),NXX,NINTMX)
        IF(NXX.EQ.0) GO TO 400
        NINT=IABS(NXX)
        IF(NINT.GT.NINTMX) THEN
          WRITE(IW,*) 'NTNFCK: CONFUSION WITH INTEGRAL FILE'
          CALL ABRT
        END IF
      ENDIF
      DO 300 M=1,NINT
C
                       NPACK = NPACKIC+M
                       IF (LABSIZ .EQ. 2) THEN
*I32                     LABEL = IX( 2*NPACK - 1 )
*I32                     IPACK = ISHFT( LABEL, -16 )
*I32                     JPACK = IAND(  LABEL, 65535 )
*I32                     LABEL = IX( 2*NPACK     )
*I32                     KPACK = ISHFT( LABEL, -16 )
*I32                     LPACK = IAND(  LABEL, 65535 )
*I64                     LABEL = IX(NPACK)
*I64                     IPACK = ISHFT( LABEL, -48 )
*I64                     JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                     KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64                     LPACK = IAND( LABEL, 65535 )
                       ELSE IF (LABSIZ .EQ. 1) THEN
*I32                     LABEL = IX(NPACK)
*I32                     IPACK = ISHFT( LABEL, -24 )
*I32                     JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32                     KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I32                     LPACK = IAND( LABEL, 255 )
*I64                     IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                       LABEL = IX( NPACK/2 )
*I64                       IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64                       JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                       LPACK = IAND( LABEL, 255 )
*I64                     ELSE
*I64                       LABEL = IX( (NPACK/2)+1 )
*I64                       IPACK = ISHFT( LABEL, -56 )
*I64                       JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64                       LPACK = IAND( ISHFT( LABEL, -32 ), 255 )
*I64                     END IF
                       END IF
                       I = IPACK
                       J = JPACK
                       K = KPACK
                       L = LPACK
C
         NIJ=IA(I)+J
         NKL=IA(K)+L
         NIK=IA(I)+K
         NIL=IA(I)+L
         NJK=IA(MAX(J,K))+MIN(J,K)
         NJL=IA(MAX(J,L))+MIN(J,L)
         VAL=XX(NPACK)
         VAL4=(VAL+VAL)+(VAL+VAL)
         F(NIJ)=F(NIJ)+VAL4*D(NKL)
         F(NKL)=F(NKL)+VAL4*D(NIJ)
         F(NIK)=F(NIK)-VAL*D(NJL)
         F(NIL)=F(NIL)-VAL*D(NJK)
         F(NJK)=F(NJK)-VAL*D(NIL)
         F(NJL)=F(NJL)-VAL*D(NIK)
  300 CONTINUE
      IF(NPACKIC.EQ.0) NPACKIC=NINTIC
C     SWITCH TO INTEGRALS ON DISK
      IF(NXX.GT.0) GO TO 100
C
  400 CONTINUE
      IF(NUM.GT.1) THEN
         DO 420 M=2,NUM
            NMAX=M-1
            DO 410 N=1,NMAX
               NIJ=IA(M)+N
               F(NIJ)=F(NIJ)*PT5
  410       CONTINUE
  420    CONTINUE
      END IF
      CALL SEQREW(IJKO)
      GO TO 2000
C
C     ----- INTEGRALS IN SQUARE ORDER -----
C
 1000 CONTINUE
      IF(NINTIC.NE.0) THEN
        WRITE(6,*) 'NTNFCK: NINTIC IS NOT READY YET'
        CALL ABRT
      ENDIF
      CALL SEQREW(IJKO)
      NA = N2
      NB = 0
      MB = 1
      ISQ= 1
      JSQ= 1
      IB = 1
C
      IJ=0
      DO 1300 I=1,NUM
      DO 1300 J=1,I
      IJ=IJ+1
      IANDJ=I.EQ.J
      CALL VCLR(A,1,NA)
C
 1040 CONTINUE
      IF (IB.LE.NB) GO TO 1060
      IF (MB.LE.0) GO TO 1100
      CALL PREAD(IJKO,XX,IX,MB,NINTMX)
      NB = IABS(MB)
      IB = 1
 1060 CONTINUE
C
      NPACK = IB
      IF (LABSIZ .EQ. 2) THEN
*I32    LABEL = IX( 2*NPACK - 1 )
*I32    IPACK = ISHFT( LABEL, -16 )
*I32    JPACK = IAND(  LABEL, 65535 )
*I32    LABEL = IX( 2*NPACK     )
*I32    KPACK = ISHFT( LABEL, -16 )
*I32    LPACK = IAND(  LABEL, 65535 )
*I64    LABEL = IX(NPACK)
*I64    IPACK = ISHFT( LABEL, -48 )
*I64    JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64    KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64    LPACK = IAND( LABEL, 65535 )
      ELSE IF (LABSIZ .EQ. 1) THEN
*I32    LABEL = IX(NPACK)
*I32    IPACK = ISHFT( LABEL, -24 )
*I32    JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32    KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I32    LPACK = IAND( LABEL, 255 )
*I64    IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64      LABEL = IX( NPACK/2 )
*I64      IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64      JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64      KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64      LPACK = IAND( LABEL, 255 )
*I64    ELSE
*I64      LABEL = IX( (NPACK/2)+1 )
*I64      IPACK = ISHFT( LABEL, -56 )
*I64      JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64      KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64      LPACK = IAND( ISHFT( LABEL, -32 ), 255 )
*I64    END IF
      END IF
      II = IPACK
      JJ = JPACK
      KK = KPACK
      LL = LPACK
C
      IF (II .NE. ISQ .OR. JJ .NE. JSQ) GO TO 1080
         KL = IA(KK) + LL
         A(KL) = XX(IB)
         IB = IB + 1
         GO TO 1040
 1080 CONTINUE
      JSQ = JSQ + 1
      IF (JSQ .LE. ISQ) GO TO 1100
      ISQ = ISQ + 1
      JSQ = 1
C
 1100 CONTINUE
      KL=0
      DO 1210 K=1,I
         LMAX=K
         IF(K.EQ.I) LMAX=J
         DO 1200 L=1,LMAX
            KL=KL+1
            KANDL=K.EQ.L
            VAL=A(KL)
            IF(VAL.EQ.ZERO) GO TO 1200
            IF(IANDJ) VAL=VAL*PT5
            IF(KANDL) VAL=VAL*PT5
            IF(IJ.EQ.KL) VAL=VAL*PT5
            VALJ=(VAL+VAL)
            VALK=-VAL*PT5
            VALJ=(VALJ+VALJ)
            VALK=(VALK+VALK)
            NIJ=IA(I)+J
            NKL=IA(K)+L
            NIK=IA(I)+K
            NIL=IA(I)+L
            NJK=IA(MAX(J,K))+MIN(J,K)
            NJL=IA(MAX(J,L))+MIN(J,L)
            F(NIJ)=F(NIJ)+VALJ*D(NKL)
            F(NKL)=F(NKL)+VALJ*D(NIJ)
            F(NIK)=F(NIK)+VALK*D(NJL)
            F(NIL)=F(NIL)+VALK*D(NJK)
            F(NJK)=F(NJK)+VALK*D(NIL)
            F(NJL)=F(NJL)+VALK*D(NIK)
 1200    CONTINUE
 1210 CONTINUE
 1300 CONTINUE
      CALL SEQREW(IJKO)
C
C     DONE CONSTRUCTING MATRIX, OFF DIAGONAL ELEMENTS NEED HALVING
C
      CALL DSCAL(N2,PT5,F,1)
      II = 0
      DO 1400 I=1,NUM
         II = II + I
         F(II) = F(II)+F(II)
 1400 CONTINUE
C
C     ----- SYMMETRIZE AND ADD BARE NUCLEUS HAMILTONIAN -----
C     IF WE HAVE A DISTRIBUTED AO INTEGRAL LIST, WE MUST GATHER F
C
 2000 CONTINUE
      IF(ITRFAO.EQ.2) CALL DDI_GSUMF(2020,F,N2)
      CALL SYMH(F,D,IA)
      CALL DAREAD(IDAF,IODA,D,N2,11,0)
      CALL VADD(F,1,D,1,F,1,N2)
      RETURN
      END
C*MODULE MCSCF   *DECK NTNONE
      SUBROUTINE NTNONE(E,H,X,G,IA,IB,IC,NE,NH,NH2,NORB,M2,NORBS,
     *                  L2,NIA,NFTX,NFTG,NFTT,ENERGY)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION E(NE),H(NH2),X(M2),G(L2),IA(NIA),IB(NORB),
     *          IC(NORB,NORBS)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00,PT5=0.5D+00)
C
C     ----- READ 1E-INTEGRALS AND -DM1- -----
C
      CALL SEQREW(NFTX)
      IF (MASWRK) THEN
         READ(NFTX) X
         READ(NFTG) G
      END IF
      CALL DAWRIT(IDAF,IODA,G,L2,68,0)
C
      DO 5 I=1,NORB
    5 IB(I)=(I-1)*NORBS
      DO 10 I=1,NIA
   10 IA(I)=(I*(I-1))/2
C
C     ----- 1E-CONTRIBUTION TO LAGRANGIAN MATRIX -----
C
      IF (MASWRK) THEN
      DO 130 IX=1,NORB
         DO 120 IG=1,NORBS
            VAL=ZERO
            DO 110 KX=1,NORBS
               I=MAX(IX,KX)
               K=MIN(IX,KX)
               IKX=IA(I)+K
               I=MAX(IG,KX)
               K=MIN(IG,KX)
               IKG=IA(I)+K
               VAL=VAL+X(IKX)*G(IKG)
  110       CONTINUE
            IIE=IB(IX)+IG
            E(IIE)=E(IIE)+VAL
  120    CONTINUE
  130 CONTINUE
      END IF
C
      IF (GOPARR) CALL DDI_GSUMF(2000,E,NE)
C
C     ----- CALCULATE ENERGY -----
C
      IF (MASWRK) THEN
      ENERGY=ZERO
      DO 210 I=1,NORBS
         IIE=IB(I)+I
         ENERGY=ENERGY+E(IIE)
         DO 200 J=1,NORBS
            II=MAX(I,J)
            JJ=MIN(I,J)
            IJ=IA(II)+JJ
            ENERGY=ENERGY+X(IJ)*G(IJ)
  200    CONTINUE
  210 CONTINUE
      ENERGY=ENERGY*PT5
C
C     ----- 1E-CONTRIBUTION TO HESSIAN MATRIX -----
C
      DO 340 IX=1,NORB
         DO 330 IG=1,NORBS
            IIH=IC(IX,IG)
            IF(IIH.EQ.0) GO TO 330
            DO 320 KX=1,IX
               KGMAX=NORBS
               IF(KX.EQ.IX) KGMAX=IG
               DO 310 KG=1,KGMAX
                  KKH=IC(KX,KG)
                  IF(KKH.EQ.0) GO TO 310
                  IH =MAX(IIH,KKH)
                  KH =MIN(IIH,KKH)
                  IKH=IA(IH)+KH
                  IKX=IA(IX)+KX
                  I=MAX(IG,KG)
                  K=MIN(IG,KG)
                  IKG=IA(I)+K
                  VAL=X(IKX)*G(IKG)
                  IF(IG.GT.IX) VAL=-VAL
                  IF(KG.GT.KX) VAL=-VAL
                  IF(KX.NE.IX.OR.KG.NE.IG) VAL=VAL+VAL
                  H(IKH)=H(IKH)+VAL
  310          CONTINUE
  320       CONTINUE
  330    CONTINUE
  340 CONTINUE
      END IF
      IF (GOPARR) THEN
         CALL DDI_BCAST(2005,'F',ENERGY,1,MASTER)
         CALL GSUMCP(2010,H,NH2)
      END IF
C
C     ----- WRITE -H- MATRIX ON TEMPORARY FILE -----
C
      CALL SEQREW(NFTT)
      IH=1
      DO 410 I=1,NH
         CALL SQWRIT (NFTT, H(IH), I)
         IH=IH+I
  410 CONTINUE
      RETURN
      END
C*MODULE MCSCF   *DECK NTNRGY
      SUBROUTINE NTNRGY(E,H,IA,IB,NE,NH,NH2,NORB,NORBS,
     *                  NIA,NFTT,ENERGY)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION E(NE),H(NH2),IA(NIA),IB(NORB)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, PT5=0.5D+00)
C
      DO 5 I=1,NORB
    5 IB(I)=(I-1)*NORBS
      DO 10 I=1,NIA
   10 IA(I)=(I*(I-1))/2
C
      IF (GOPARR) CALL DDI_GSUMF(2000,E,NE)
C
C     ----- CALCULATE ENERGY -----
C
      IF (MASWRK) THEN
         ENERGY=ZERO
         DO 210 I=1,NORBS
            IIE=IB(I)+I
            ENERGY=ENERGY+E(IIE)
  210    CONTINUE
         ENERGY=ENERGY*PT5
      END IF
      IF (GOPARR) THEN
         CALL DDI_BCAST(2005,'F',ENERGY,1,MASTER)
         CALL GSUMCP(2010,H,NH2)
      END IF
C
C     ----- WRITE -H- MATRIX ON TEMPORARY FILE -----
C
      CALL SEQREW(NFTT)
      IH=1
      DO 410 I=1,NH
         CALL SQWRIT (NFTT, H(IH), I)
         IH=IH+I
  410 CONTINUE
      RETURN
      END
C*MODULE MCSCF   *DECK NTNOUT
      SUBROUTINE NTNOUT(E,M,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION E(M,1)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
 9999 FORMAT(/)
 9998 FORMAT(I5,8F15.10)
      MAX=8
      IMAX=0
  100 IMIN=IMAX+1
      IMAX=IMAX+MAX
      IF(IMAX.GT.M) IMAX=M
      IF (MASWRK) WRITE(IW,9999)
      DO 400 J=1,N
      IF (MASWRK) WRITE(IW,9998) J,(E(I,J),I=IMIN,IMAX)
  400 CONTINUE
      IF(IMAX.LT.M) GO TO 100
      RETURN
      END
C*MODULE MCSCF   *DECK NTNPAR
      SUBROUTINE NTNPAR(NROT,NROT2,IB,IC,LABMO,NOROT,MXNORO,MOFRZ,
     *                  NORB,NORBS,NCORBS,FORS,FCORE,DBUG)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL FORS,FCORE,CORCOR,ACTACT,ACTCOR,VIRCOR
      LOGICAL DBUG
C
      DIMENSION IB(NORB),IC(NORB,NORBS),LABMO(NORBS),NOROT(2,MXNORO),
     *          MOFRZ(15)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
C
C     ----- INITIALIZE ORBITAL ROTATION PARAMETER ARRAY -----
C     ZERO ELEMENTS OF -IC- INDICATE REDUNDANT ROTATIONS, OR
C     VARIATIONAL PARAMETERS WHICH ARE NOT TO BE OPTIMIZED.
C     -NOROT- ALLOWS FOR USER INPUT OF IGNORABLE ROTATIONS.
C
      DO 100 I=1,NORB
         IB(I)=(I-1)*NORBS
  100 CONTINUE
C
C        USE OF SYMMETRY TO ELIMINATE ROTATIONS CAN BE VERY EFFECTIVE.
C        FOR THIS WE NEED THE MO IRRREP SYMMETRY LABELS FROM -TRFSYM-
C
      CALL DAREAD(IDAF,IODA,LABMO,NORB,324,1)
C
      NROT=0
      DO 230 I=1,NORB
         LABI = LABMO(I)
         DO 202 IFRZ=1,15
            IF(I.EQ.MOFRZ(IFRZ)) GO TO 230
  202    CONTINUE
         JMAX=MIN(I,NORBS)
         DO 220 J=1,JMAX
            IC(I,J)=0
            IF(I.LE.NORBS) IC(J,I)=0
C
            IF(J.EQ.I)                          GO TO 220
C
            CORCOR=(I.LE.NCORBS)                .AND.
     *             (J.LE.NCORBS)
            IF(CORCOR)                          GO TO 220
C
            ACTACT=(I.GT.NCORBS.AND.I.LE.NORBS) .AND.
     *             (J.GT.NCORBS.AND.J.LE.NORBS)
            IF(FORS.AND.ACTACT)                 GO TO 220
C
C              ORBITALS OF DIFFERENT SYMMETRY SHOULD NOT MIX
C
            IF(LABI.NE.LABMO(J))                GO TO 220
C
            ACTCOR=(I.GT.NCORBS.AND.I.LE.NORBS) .AND.
     *             (J.LE.NCORBS)
            VIRCOR=(I.GT.NORBS)                 .AND.
     *             (J.LE.NCORBS)
            IF(FCORE.AND.(ACTCOR.OR.VIRCOR))    GO TO 220
C
C              NEXT TWO LOOPS CHECK FOR ROTATIONS OR MOS FROZEN BY USER
C
            DO 210 INOROT=1,MXNORO
               II=NOROT(1,INOROT)
               JJ=NOROT(2,INOROT)
               IF(II.EQ.I.AND.JJ.EQ.J) GO TO 220
  210       CONTINUE
C
            DO 212 IFRZ=1,15
               IF(J.EQ.MOFRZ(IFRZ)) GO TO 220
  212       CONTINUE
C
C             FOUND AN ORBITAL ROTATION THAT MUST BE OPTIMIZED
C
            NROT=NROT+1
            IC(I,J)=NROT
            IF(I.LE.NORBS) IC(J,I)=NROT
C
  220    CONTINUE
  230 CONTINUE
      NROT2=(NROT*(NROT+1))/2
C
      IF(.NOT.DBUG) RETURN
      WRITE(IW,9999)
      DO 300 I=1,NORB
         WRITE(IW,9998) (IC(I,J),J=1,NORBS)
  300 CONTINUE
      RETURN
C
 9999 FORMAT(1X,'ORBITAL PAIR ROTATION PARAMETERS = ')
 9998 FORMAT(10X,40I3)
      END
C*MODULE MCSCF   *DECK NTNROT
      SUBROUTINE NTNROT(DEMAX,IROTMX,JROTMX,E,IROT,IC,NROT,NORB,NORBS,
     *                  NSKIP,ACURCY,SOME,OUT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SOME,OUT
C
      DIMENSION E(NORBS,NORB),IROT(NORBS,NORB),IC(NORB,NORBS)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- GET LAGRANGIAN ASYMMETRY -----
C     DEMAX IS MAX ASYM, -IC- HOLDS ROTATIONS BEING OPTIMIZED.
C     THE TRADITIONAL DEFAULT FOR -TOL- WAS 1D-8, FOR MANY YEARS.
C
      TOL = MAX(1.0D-09,1.0D-03*ACURCY)
      DO 100 I=1,NORB
      DO 100 J=1,NORBS
  100 IROT(J,I)=0
      NROT=0
      DEMAX=ZERO
      DO 220 I=2+NSKIP,NORB
         JMAX=MIN(I-1,NORBS)
         DO 210 J=1+NSKIP,JMAX
            DUM=E(J,I)
            IF(I.LE.NORBS) DUM=DUM-E(I,J)
            DUM= ABS(DUM)
            IF(DUM.GT.DEMAX) THEN
               IROTMX=I
               JROTMX=J
               DEMAX=DUM
            END IF
            IROT(J,I)=0
            IF(I.LE.NORBS) IROT(I,J)=0
            IF(DUM.LT.TOL) GO TO 210
            IJ=IC(I,J)
            IF(IJ.EQ.0) THEN
               IF(SOME) WRITE(IW,9996) DUM,I,J
            ELSE
               NROT=NROT+1
               IROT(J,I)=NROT
               IF(I.LE.NORBS) IROT(I,J)=NROT
            END IF
  210    CONTINUE
  220 CONTINUE
      IF(SOME) WRITE(IW,9997) NROT
C
      IF(.NOT.OUT) RETURN
      WRITE(IW,9999)
      DO 300 I=1,NORB
         WRITE(IW,9998) (IROT(J,I),J=1,NORBS)
  300 CONTINUE
      RETURN
C
 9999 FORMAT(' ORBITAL PAIR ROTATION PARAMETERS = ')
 9998 FORMAT(10X,40I3)
 9997 FORMAT(' NUMBER OF INDEPENDENT ORBITAL ROTATION PARAMETERS = ',I5)
 9996 FORMAT(1X,'GRAD=',1P,E14.6,0P,' FOR ORB. ROTATION I,J=',
     *          2I5,' EXCEEDS TOLERANCE.')
      END
C
C*MODULE MCSCF   *DECK NTNRPH
      SUBROUTINE NTNRPH(IROT,IC,E,A,IA,H,B, AD,V,W,T
     *,                 A2,AA,VEC,EIG,WRK,IWRK
     *,                 NROT,NROT1,NROT2,NIA,NORB,NORBS
     *,                 NFTT,OUT,SQCDF,ACURCY,MXXPAN)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL OUT,DBUG
C
      DIMENSION IROT(NORBS,NORB),IC(NORB,NORBS),E(NORBS,NORB)
     *,         A(NROT2),IA(NIA),H(*),B(NROT)
     *,         AD(NROT1),V(NROT1,MXXPAN),W(NROT1,MXXPAN),T(NROT1)
     *,         A2(*),AA(*),VEC(MXXPAN,MXXPAN),EIG(MXXPAN)
     *,         WRK(MXXPAN,8),IWRK(MXXPAN)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
C
      PARAMETER (ZERO=0.0D+00, PT5=0.5D+00, HSHIFT=0.5D+00,
     *           TEST=0.5D+00, TEST1=0.1D+00)
C
      DO 5 I=1,NIA
    5 IA(I)=(I*(I-1))/2
      CALL VCLR(A,1,NROT2)
C
C     ----- LAGRANGIAN CONTRIBUTION TO -A- -----
C
      DO 140 IX=1,NORB
         DO 130 IG=1,NORBS
            IB=IROT(IG,IX)
            IF(IB.EQ.0) GO TO 110
               VAL=E(IG,IX)
               IF(IG.GT.IX) VAL=-VAL
               IJ=IA(IB+1)+1
               A(IJ)=A(IJ)+(VAL+VAL)
  110       CONTINUE
C
            DO 120 K=1,NORB
               IF(IX.GT.NORBS.AND.K.GT.NORBS) GO TO 120
               VAL=E(IG,IX)
               IF(K.GT.IX) VAL=-VAL
               IF(K.LT.IG) VAL=-VAL
               IF(IX.GT.NORBS) IB=IROT(K,IX)
               IF(IX.LE.NORBS) IB=IROT(IX,K)
               JB=IROT(IG,K)
               IF(IB.EQ.0.OR.JB.EQ.0) GO TO 120
                  I=MAX(IB,JB)
                  J=MIN(IB,JB)
                  IJ=IA(I+1)+(J+1)
                  A(IJ)=A(IJ)+VAL
  120       CONTINUE
  130    CONTINUE
  140 CONTINUE
C
C     ----- OTHER CONTRIBUTIONS TO -A- -----
C
      CALL SEQREW(NFTT)
      DO 250 IX=2,NORB
         MAXIG=MIN(IX-1,NORBS)
         DO 240 IG=1,MAXIG
            IIC=IC(IX,IG)
            IF(IIC.EQ.0) GO TO 240
            CALL SQREAD (NFTT, H(1), IIC)
            IJH=0
            IIA=IROT(IG,IX)
            DO 220 JX=2,IX
               MAXJG=MIN(JX-1,NORBS)
               IF(JX.EQ.IX) MAXJG=IG
               DO 210 JG=1,MAXJG
                  JJC=IC(JX,JG)
                  IF(JJC.EQ.0) GO TO 210
                  IJH=IJH+1
                  JJA=IROT(JG,JX)
                  IF(IIA.EQ.0.OR.JJA.EQ.0) GO TO 210
                  IJA=IA(IIA+1)+(JJA+1)
                  A(IJA)=A(IJA)+H(IJH)
  210          CONTINUE
  220       CONTINUE
  240    CONTINUE
  250 CONTINUE
C
C     ----- MULTIPLY DIAGONAL ELEMENTS BY 2 -----
C
      DO 270 I=1,NROT
         II=IA(I+1)+(I+1)
         A(II)=A(II)+A(II)
  270 CONTINUE
C
C     ----- PRINT AUGMENTED HESSIAN MATRIX -----
C
      IF(OUT) THEN
         WRITE(IW,9999)
         DO 280 I=1,NROT
            N=IA(I+1)+1
            WRITE(IW,9998) (A(N+J),J=1,I)
            WRITE(IW,9997)  A(N)
  280    CONTINUE
      END IF
C
C     ----- B. LENGSFIELD'S LEVEL SHIFT -----
C
      DO 290 I=1,NROT
         N=IA(I+1)+(I+1)
         IF(A(N).LE.ZERO) A(N)=HSHIFT
         AD(I+1)=A(N)
  290 CONTINUE
      AD(1) = ZERO
C
C     ----- FIND LOWEST ROOT OF AUGMENTED HESSIAN MATRIX -----
C     AN IN-MEMORY DAVIDSON SOLVER IS INVOKED TO GET THE LOWEST ROOT
C
      MAXIT  = 200
      CVGTOL = ACURCY
      PRTTOL = 1.0D-09
      DBUG   = .FALSE.
      IDUM=1
      DUM=0.0D+00
      CALL NTNDVD(A,AD,V,W,T, AA,A2,VEC,EIG,WRK,IWRK, IA,
     *            NROT1,MXXPAN,CVGTOL,MAXIT,OUT,DBUG,PRTTOL,
     *            IDUM,IDUM,IDUM,IDUM,DUM,DUM,DUM,DUM,DUM,IDUM,DUM)
C
C     WE COULD CHECK THE DAVIDSON CODE WITH GLDIAG, IF DESIRED
C
C-----CALL GLDIAG(NROT1,1,NROT1,A,W,AD,V,IERR,T)
C
      IF(OUT) WRITE(IW,9996) T(1),(V(I,1),I=1,NROT1)
C
C     ----- B. LENGSFIELD'S EMPIRICAL FUDGING -----
C
      TESTV =TEST
      TESTV1=TEST1
      TESTV2= SQRT(TEST)
      C1=V(1,1)
      IF( ABS(C1).GT.TESTV2) GO TO 340
      C1SQ=C1**2
C
  300 CONTINUE
      C2SQ=ZERO
      DO 310 I=1,NROT
         C2=V(I+1,1)
         IF( ABS(C2).LT.TESTV1) GO TO 310
         C2SQ=C2SQ+C2**2
  310 CONTINUE
      IF(OUT) WRITE(IW,9994) C1SQ,C2SQ,TESTV1
      IF(C1SQ+C2SQ.GT.TESTV) GO TO 320
      TESTV1=TESTV1*PT5
      IF(OUT) WRITE(IW,9993)
      GO TO 300
C
  320 CONTINUE
      V(1,1)= SIGN(TESTV2,V(1,1))
      SCALE= SQRT((C1SQ+C2SQ-TESTV)/C2SQ)
      DO 330 I=1,NROT
         C2=V(I+1,1)
         IF( ABS(C2).LT.TESTV1) GO TO 330
         V(I+1,1)=C2*SCALE
  330 CONTINUE
C
  340 CONTINUE
      IF(V(1,1).GT.ZERO) GO TO 360
      DO 350 I=1,NROT1
         V(I,1)=-V(I,1)
  350 CONTINUE
C
  360 CONTINUE
      SQCDF=ZERO
      DO 370 I=1,NROT
         DUM=V(I+1,1)
         SQCDF=SQCDF+DUM*DUM
         B(I)=DUM
  370 CONTINUE
      IF(OUT) WRITE(IW,9995) SQCDF
C
C     ----- CONSTRUCT ROTATION PARAMETERS IN -E- -----
C
      DO 420 IX=1,NORB
         DO 410 IG=1,NORBS
            E(IG,IX)=ZERO
            IB=IROT(IG,IX)
            IF(IB.EQ.0) GO TO 410
            VAL=B(IB)
            IF(IG.GT.IX) VAL=-VAL
            E(IG,IX)=VAL
  410    CONTINUE
  420 CONTINUE
      RETURN
C
 9999 FORMAT(40X,'   HESSIAN MATRIX   ',40X,' GRADIENT ',/)
 9998 FORMAT(5E20.12)
 9997 FORMAT(100X,E20.12,/)
 9996 FORMAT(' EIGENVECTOR OF AUGMENTED HESSIAN ',F20.12,/,(5F20.12))
 9995 FORMAT(' SQCDF = ',E20.4)
 9994 FORMAT(' VECTOR FIXUP USED. C1SQ = ',F10.8,' C2SQ = ',F10.8,
     1 ' TESTV1 = ',F7.3)
 9993 FORMAT(' TESTV1 REDUCED BY FACTOR OF 2. ')
      END
C*MODULE MCSCF   *DECK NTNTRF
      SUBROUTINE NTNTRF(NUM,NORB,NORBS,NCORBS,FCORE,DAMP,S,V,U,T)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL FIRST,FCORE
C
      DIMENSION S(NORBS,NORB),V(NUM,NUM),U(NORB,NORB),T(NORB)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00,
     *           ROTMAX=0.2D+00)
C
C     ----- APPROXIMATE ORBITAL TRANSFORMATION -----
C
C          S      = ROTATION PARAMETERS
C          U      = EXP(S) EXPANDED TO SECOND ORDER
C          U      = I + S * ( I + S/2 ) + ...
C
C     NOTE - IF ANY OF THE S(I,J) IS GREATER THAN 0.2 IN ABSOLUTE
C     VALUE, EXPAND EXP(S) TO FIRST ORDER ONLY.
C
      FIRST=.FALSE.
      DO 120 I=2,NORB
         JMAX=MIN(I-1,NORBS)
         DO 110 J=1,JMAX
            IF( ABS(S(J,I)).LT.ROTMAX) GO TO 110
            FIRST=.TRUE.
  110    CONTINUE
  120 CONTINUE
C
      CALL VCLR(U,1,NORB*NORB)
      DO 210 I=1,NORB
         U(I,I)=ONE
  210 CONTINUE
      IF(FIRST) GO TO 300
C
      FAC=ONE/TWO
      FAC=FAC/(ONE+DAMP)
      DO 250 I=2,NORB
         JMAX=MIN(I-1,NORBS)
         DO 240 J=1,JMAX
            DUM=S(J,I)*FAC
            U(I,J)= DUM
            U(J,I)=-DUM
  240    CONTINUE
  250 CONTINUE
C
  300 CONTINUE
      DO 320 I=1,NORB
         DO 310 J=1,NORB
            V(I,J)=ZERO
  310    CONTINUE
  320 CONTINUE
C
      FAC=ONE/(ONE+DAMP)
      DO 360 I=2,NORB
         JMAX=MIN(I-1,NORBS)
         DO 350 J=1,JMAX
            DUM=S(J,I)*FAC
            DO 340 K=1,NORB
               V(J,K)=V(J,K)-DUM*U(I,K)
               V(I,K)=V(I,K)+DUM*U(J,K)
  340       CONTINUE
  350    CONTINUE
  360 CONTINUE
C
      DO 390 I=1,NORB
         DO 380 J=1,NORB
            U(I,J)=V(I,J)
            V(I,J)=ZERO
  380    CONTINUE
         U(I,I)=U(I,I)+ONE
  390 CONTINUE
C
C     ----- OPTION TO FREEZE CORE ORBITALS -----
C     THIS IS TO MAKE DOUBLY CERTAIN U DOESN'T ROTATE CORE ORBS,
C     THESE PARTS OF THE MATRIX ARE PRBABLY ZERO/ONE ALREADY.
C
      IF(FCORE) THEN
         DO 420 I=1,NCORBS
            DO 410 J=1,NORB
               U(I,J)=ZERO
               U(J,I)=ZERO
  410       CONTINUE
            U(I,I)=ONE
  420    CONTINUE
      END IF
C
C     ----- ORTHONORMALIZE THE TRANSFORMATION MATRIX -----
C
      DO 560 I=1,NORB
         DUM=ZERO
         DO 510 K=1,NORB
            DUM=DUM+U(K,I)*U(K,I)
  510    CONTINUE
         DUM= ONE/ SQRT(DUM)
         DO 520 K=1,NORB
            U(K,I)=U(K,I)*DUM
  520    CONTINUE
         IF(I.EQ.NORB) GO TO 560
         I1=I+1
         DO 550 J=I1,NORB
            DUM=ZERO
            DO 530 K=1,NORB
               DUM=DUM+U(K,I)*U(K,J)
  530       CONTINUE
            DO 540 K=1,NORB
               U(K,J)=U(K,J)-DUM*U(K,I)
  540       CONTINUE
  550    CONTINUE
  560 CONTINUE
C
C     ----- ROTATE THE ORBITALS -V- BY TRANSFORMATION -U- -----
C
      N3=NUM*NUM
      CALL DAREAD(IDAF,IODA,V,N3,15,0)
      DO 650 I=1,NUM
         DO 630 J=1,NORB
            DUM=ZERO
            DO 620 K=1,NORB
               DUM=DUM+V(I,K)*U(K,J)
  620       CONTINUE
            T(J)=DUM
  630    CONTINUE
         DO 640 J=1,NORB
            V(I,J)=T(J)
  640    CONTINUE
  650 CONTINUE
      CALL DAWRIT(IDAF,IODA,V,N3,15,0)
      RETURN
      END
C*MODULE MCSCF   *DECK NTNDE1
      SUBROUTINE NTNDE1(V,DCOR,DVAL,GIJ,IA,NORBS,NCORBS,HO,E,NIA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      DIMENSION GIJ(*),V(NUM,*),IA(*),DCOR(*),DVAL(*),HO(*)
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00)
C
      DO 5 I=1,NIA
         IA(I)=(I*I-I)/2
    5 CONTINUE
C
C     ----- FORM TOTAL DENSITY MATRIX IN -AO- BASIS -----
C
C     M AND N ARE AO INDICES, K AND L ARE MO INDICES
C
C     DCOR(M,N) = SUM 2 * CMK * CNK OVER K CORE ORBITALS
C
C     DVAL(M,N) = SUM GKL * CMK * CNL OVER K,L ACTIVE ORBITALS
C
C     ECORE = SUM (DCORE(M,N) + DVAL(M,N)) * H(M,N) OVER M,N
C
C
      E = ZERO
      N2=NUM*(NUM+1)/2
      NCI= NORBS-NCORBS
      N3=NUM*NUM
      CALL DAREAD(IDAF,IODA,V,N3,15,0)
C
C     ----- CONSTRUCTION DCOR -----
C
      MN=0
      DO 10 M=1,NUM
      DO 10 N=1,M
      MN=MN+1
   10 DCOR(MN)=ZERO
C
      IF(NCORBS.EQ.0) GO TO 40
C
      MN=0
      DO 30 M=1,NUM
      DO 30 N=1,M
      DUM=ZERO
      DO 20 I=1,NCORBS
   20 DUM=DUM + TWO*V(M,I)*V(N,I)
      MN=MN+1
      DCOR(MN)=DUM
   30 CONTINUE
C
   40 CONTINUE
C
C     ----- CONSTRUCTION DVAL -----
C
      MN=0
      DO 70 M=1,NUM
      DO 70 N=1,M
      DUM = ZERO
      DO 60 K= 1,NCI
      DO 50 L= 1,K
      KL=IA(K)+ L
      GMN =  GIJ(KL) * V(M,K+NCORBS) * V(N,L+NCORBS)
      GNM =  GIJ(KL) * V(N,K+NCORBS) * V(M,L+NCORBS)
   50 DUM = DUM + GMN + GNM
      DUM = DUM - GMN
   60 CONTINUE
      MN=MN+1
      DVAL(MN)=DUM
   70 CONTINUE
C
      CALL DAREAD(IDAF,IODA,HO,N2,11,0)
      DO 90 I = 1,NUM
      DO 80 J = 1,I
      IJ =IA(I) + J
      DT=DCOR(IJ)+DVAL(IJ)
   80 E = E + DT*HO(IJ)*TWO
   90 E = E - DT*HO(IJ)
C
      RETURN
      END
C*MODULE MCSCF   *DECK NTNLAG
      SUBROUTINE NTNLAG(IA,FC,FA,GIJ,EG,NORB,NORBS,NCORBS,IB,NIA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     -----  PARTIAL CONTRIBUTION TO E(RI)  -----
C
C     I CORE  EG(RI) = 2* ( FC(RI) + FA(RI))
C
C     I VAL   EG(RI) = SUM GIJ * FC(RJ) OVER J
C
      DIMENSION EG(1),GIJ(1),IA(1),IB(1)
      DIMENSION FA(1),FC(1)
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00)
C
      NCI = NORBS-NCORBS
C
      DO 5 I=1,NIA
         IA(I)=(I*I-I)/2
    5 CONTINUE
C
       IF(NCORBS.EQ.0) GO TO 30
C
C     ----- CORE CONTRIBUTION -----
C
       DO 20 M = 1,NORB
       DO 20 I = 1,NCORBS
       IM = IA(MAX0(M,I))+MIN0(M,I)
       CNT = (FC(IM) + FA(IM))*TWO
       IR1= IB(M)+ I
       EG(IR1)= CNT
   20  CONTINUE
C
C      ----- VALENCE CONTRIBUTION -----
C
   30  CONTINUE
C
       DO 60 M=1,NORB
       DO 50 I=1,NCI
       DUM=ZERO
       DO 40 J=1,NCI
       IJ=IA(MAX0(I,J))+MIN0(I,J)
       XGIJ=GIJ(IJ)
       MJ=IA(MAX0(M,J+NCORBS))+MIN0(M,J+NCORBS)
       DUM =  DUM + FC(MJ)* XGIJ
   40  CONTINUE
       IIR= IB(M)+ I + NCORBS
   50  EG(IIR)=DUM
   60  CONTINUE
C
      RETURN
      END
C*MODULE MCSCF   *DECK NTNHES
      SUBROUTINE NTNHES(IA,IC,FC,FA,GIJ,EH,
     *                  NORB,NORBS,NCORBS,NIA)
C
C     ----- PARTIAL CONTRIBUTION TO EH(RI,SJ ) -----
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION EH(1),GIJ(1),IA(1),IC(NORB,1)
      DIMENSION FC(1),FA(1)
C
      PARAMETER (PT5=0.5D+00)
C
      DO 5 I=1,NIA
         IA(I)=(I*I-I)/2
    5 CONTINUE
C
      IF (NCORBS.EQ.0) GO TO 25
C
C     CONTRIBUTION TO H(RI,JS) FROM CORE I=J
C
C     H(RI,SJ) = GIJ* ( FC(RS) + FA (RS) )
C
       DO 20 M = NCORBS+1,NORB
       DO 20 N = NCORBS+1,M
       MN=IA(M)+N
C
       DO 15 I=1,NCORBS
       NR1=IC(M,I)
       IF(NR1.EQ.0) GO TO 15
       NR2=IC(N,I)
       IF(NR2.EQ.0) GO TO 15
       MNR1= MAX0(NR1,NR2)
       MNR2= MIN0(NR1,NR2)
       MINI=IA(MNR1)+MNR2
       FF= FC(MN)+FA(MN)
       FF = FF + FF
       IF(M.NE.N) FF= FF + FF
       EH(MINI)=EH(MINI)+ FF
   15  CONTINUE
C
   20  CONTINUE
C
   25  CONTINUE
C
C     CONTRIBUTION TO H(RI,JS) FROM VALENCE I,J
C
C     H(RI,SJ) = GIJ *  FC(RS)
C
       DO 60 M = 1,NORB
       DO 60 N = 1,M
       MN=IA(M)+N
C
       DO 50 I=NCORBS+1,NORBS
        MI1=IC(M,I)
        NI1=IC(N,I)
C
         DO 40 J=NCORBS+1,I
          NJ2=IC(N,J)
          MJ2=IC(M,J)
C
          II=I-NCORBS
          JJ=J-NCORBS
          IJ=IA(II) + JJ
          XX=GIJ(IJ)*FC(MN)
          XX= XX+XX
          IF(M.EQ.N)XX= XX*PT5
          IF(I.EQ.J)XX= XX*PT5
C
          IF(MI1.EQ.0.OR.NJ2.EQ.0) GO TO 30
          FF = XX
          IF(M.LT.I)FF=-FF
          IF(N.LT.J)FF=-FF
          MINJ=IA(MAX0(MI1,NJ2))+MIN0(MI1,NJ2)
          EH(MINJ)=EH(MINJ) + FF
C
   30     IF(NI1.EQ.0.OR.MJ2.EQ.0) GO TO 40
          FF = XX
          IF(N.LT.I)FF=-FF
          IF(M.LT.J)FF=-FF
          NIMJ=IA(MAX0(NI1,MJ2))+MIN0(NI1,MJ2)
          EH(NIMJ)=EH(NIMJ) + FF
   40     CONTINUE
C
   50    CONTINUE
C
   60  CONTINUE
C
      RETURN
      END
C*MODULE MCSCF   *DECK MCAODBF
      SUBROUTINE MCAODBF(V,GMN,GIJ,NUM,NORBS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      DIMENSION GIJ(*),GMN(*),V(NUM,NUM)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- FORM TOTAL DENSITY MATRIX IN -AO- BASIS -----
C
      N3=NUM*NUM
      L2=(NORBS*NORBS+NORBS)/2
      CALL DAREAD(IDAF,IODA,GIJ,L2,68,0)
      CALL DAREAD(IDAF,IODA,V,N3,15,0)
      MN=0
      DO 740 M=1,NUM
         DO 730 N=1,M
            DUM=ZERO
            DO 720 I=1,NORBS
               DO 710 J=1,NORBS
                  II=MAX(I,J)
                  JJ=MIN(I,J)
                  IJ=IA(II)+JJ
                  DNORB=GIJ(IJ)
                  DUM=DUM+DNORB*V(M,I)*V(N,J)
  710          CONTINUE
  720       CONTINUE
            MN=MN+1
            GMN(MN)=DUM
  730    CONTINUE
  740 CONTINUE
      RETURN
      END
