C 11 DEC 06 - MWS - ADJUST COMMON BLOCKS FOR FMO CHANGES
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 25 OCT 06 - SPW - NEW MODULE FOR DETERMINANT BASED NEO-FULL CI
C
C*MODULE NEOFCI  *DECK NEOFCI
      SUBROUTINE NEOFCI(NEOCIP,NPASS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DOUBLE PRECISION NEODET,MCSCF
      DOUBLE PRECISION NEOSCF,NEOCI
C
      DIMENSION NEOCIP(10)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,DOEXCH,
     *        DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,DOVVVO,DOVVVV,DOCORE
      LOGICAL DIRSCF,FDIFF
C
      PARAMETER (MXRT=100, MXSH=5000, MXATM=2000)
C
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,S,SZ,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORB,
     *                NA,NB,K,KST,IROOT,IPURES,MAXW1,NITER,MAXP,NCI,
     *                IGPDET,KSTSYM,NFTGCI
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /NEOMTD/ NEOSCF,NEOCI
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     &                TSYM(432),INVT(48),NT
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA CIDET/8HCIDET   /, RNONE/8HNONE    /
      DATA NEODET/8HNEODET  /
      DATA MCSCF/8HMCSCF   /
C
C     --- MAIN  DRIVER FOR NEOFCI CALCULATION ---
C
      IF(NEOSCF.NE.MCSCF .AND. NPASS.EQ.0) THEN
C
C        --- IF NEO-CI BASED ON NEO-HF, ALL SYMMETRY SWITCHED
C            OFF HERE FOR NOW (EARLIER FOR MCSCF) ---
C
         IF(NT.NE.1) THEN
            IF(MASWRK) WRITE(IW,8000)
            IF(MASWRK) WRITE(IW,8010)
            NT = 1
            IF(.NOT.DIRSCF) CALL JANDK
            CALL NEOSET(.FALSE.,.FALSE.)
         END IF
C
      END IF
C
      IF(NPASS.EQ.0 .AND. MASWRK) WRITE(IW,9000)
C
C     --- DOING NEO-CI BASED ON NEO-HF ORBITALS ---
C
      IF(NEOSCF.NE.MCSCF) THEN
C
C        --- GET ELECTRONIC DETERMINANT CI INPUT ---
C
         CALL EDETIN(NEOCIP(1),CIDET)
C
C        --- GET NUCLEAR DETERMINANT CI INPUT PLUS NEO-CI OPTIONS ---
C
         CALL NDETIN(NEOCIP(1),NEODET)
C
C        --- NEED ELECTRONIC AO TO MO TRANSFORMATION ----
C
         DDITRF=GOPARR
         DOOOOO=.TRUE.
         DOVOOO=.FALSE.
         DOVVOO=.FALSE.
         DOVOVO=.FALSE.
         DOVVVO=.FALSE.
         DOVVVV=.FALSE.
         DOCORE=.TRUE.
         DOEXCH=SCFTYP.EQ.RNONE
         CALL TRFMCX(NEOCIP(2),NCOR,NORB,NORB,.FALSE.,DOEXCH,
     *               DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *               DOVVVO,DOVVVV,DOCORE)
C
      END IF
C
C     --- CARRY OUT NUCLEAR-NUCLEAR AND NUCLEAR-ELECTRONIC AO TO MO
C         TRANSFORMATIONS ---
C
      CALL NEOTRN(NPASS)
C
C     --- CARRY OUT NEO-CI CALCULATION ---
C
      CALL NUELCI(NPASS)
C
      RETURN
C
 8000 FORMAT(/1X,'**** SWITCHING OFF ALL SYMMETRY FOR NEO-CI ****'/)
 8010 FORMAT(6X,'... MAY NEED TO RECALCULATE 2-PARTICLE AO INTEGRALS')
 9000 FORMAT(/5X,54(1H-)/
     * 5X,'NUCLEAR ELECTRONIC ORBITAL (NEO) DETERMINANTAL FULL CI'/
     * 5X,'    PROGRAM WRITTEN BY SIMON P. WEBB, MIKE PAK,'/
     * 5X,'                    AND CHET SWALINA'/
     * 5X,54(1H-)///
     * 5X,'... PARTS OF THIS PROGRAM USE SUBROUTINES FROM THE'/
     * 5X,'    AMES LABORATORY DETERMINANTAL FULL CI PROGRAM'/
     * 5X,'    WRITTEN BY JOE IVANIC AND KLAUS RUEDENBERG')
C
      END
C*MODULE NEOFCI   *DECK NUELCI
      SUBROUTINE NUELCI(NPASS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DOUBLE PRECISION METHOD,NEOSCF,NEOCI,MCSCF
C
      PARAMETER (MXATM=2000)
      PARAMETER (MXNORO=250, MXRT=100, MXAO=8192)
C
      LOGICAL SOME1,SOME2
      LOGICAL GOPARR,DSKWRK,MASWRK,PACK2E
      LOGICAL CANONC,FCORE,FORS,NOCI,EKT,LINSER
      LOGICAL DOTWON,CLABEL
C
      DIMENSION TIMSTR(3)
C
      COMMON /CIFILS/ NFT11,NFT12,NFT13,NFT14,NFT15,NFT16,IDAF20,NEMMX
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,S,SZ,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORB,
     *                NA,NB,K,KST,IROOT,IPURES,MAXW1,NITER,MAXP,NCI,
     *                IGPDET,KSTSYM,NFTGCI
      COMMON /DETNEO/ WSTATE_NEO(MXRT),SPINS_NU(MXRT),CRIT_NU,
     *                PRTTOL_NEO,S_NU,SZ_NU,GRPDET_NU,STSYM_NU,
     *                NFLGDM_NEO(MXRT),IWTS_NU(MXRT),NCORSV_NU,
     *                NCOR_NU,NACT_NU,NORB_NU,NA_NU,NB_NU,
     *                K_NEO,KST_NEO,IROOT_NEO,IPURES_NEO,MAXW1_NEO,
     *                NITER_NEO,MAXP_NEO,NCI_NU,IGPDET_NU,KSTSYM_NU
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZE,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FUNCT / E,EGRAD(3,MXATM)
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NAINF,NBINF,
     *                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 /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,
     *                NORBMC,NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /NEOCIR/ IEROOT
      COMMON /NEOMTD/ NEOSCF,NEOCI
      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      COMMON /OUTPUT/ NPRINTX,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA ALDET/8HALDET   /
      DATA MCSCF/8HMCSCF   /
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
C
C     --- DRIVER TO DO MIXED NUCLEAR/ELECTRONIC MOLECULAR ORBITAL (NEO)
C         FULL CONFIGURATION INTERACTION (FCI) WITHIN GIVEN ELECTRONIC
C         AND NUCLEAR ACTIVE SPACES ---
C
C         PLUS: OPTIONAL NEO-CI DENSITY MATRICES (ELECTRONIC AND NUCLEAR
C               1- AND 2-PARTICLE, AS WELL AS MIXED NUC/ELEC).
C         PLUS: OPTIONAL NEO NATURAL ORBITALS ELECTRONIC AND NUCLEAR.
C
C     --- NEO-FCI CODE WRITTEN BY SIMON P. WEBB, MIKE PAK, AND
C         CHET SWALINA ---
C         (IT IS MAINLY A CLEANED UP VERSION OF CODE WRITTEN IN 2001-2002
C          BY SIMON P. WEBB)
C
      SOME1 = MASWRK .AND. NPASS.EQ.0
      SOME2 = MASWRK .AND.  (NPASS.EQ.0 .OR. NPASS.EQ.-1)
C
      CLABEL = .FALSE.
C
C     --- SET IROOT IF MCSCF STATE AVERAGING WITH
C              MCSCF OUTPUT ENERGY A PARTICULAR STATE ---
C
      IF(IEROOT.GT.0) IROOT = IEROOT
C
C     --- NEO MUST USE DETERMINANTS ---
C
      CISTEP = ALDET
C
C     --- INITIALIZE SOME ELECTRONIC VARIABLES ---
C
      NSYM = 2**IGPDET
C
C     - ELECTRONIC BASIS FUNCTIONS
      L0 = NQMT
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
C
C     - ELECTRONIC ACTIVE SPACE
      NTOT = NACT + NCORSV       ! ACTIVE ORBS + CORE ORBS
      NTOT2 = (NTOT*NTOT+NTOT)/2
      NORB = NACT
C
C     --- NOW INITIALIZE SOME NUCLEAR CI VARIABLES ---
C
C     - NUCLEAR BASIS FUNCTIONS
      N1 = NUMNB
      N2 = (N1*N1+N1)/2
      N3 = N1*N1
C
C     - NUCLEAR ACTIVE SPACE
      NTOT_NU = NACT_NU + NCORSV_NU           ! ACTIVE NUC ORBS + CORES
      NTOT2_NU = (NTOT_NU*NTOT_NU+NTOT_NU)/2
      NORB = NACT
      NORB_NU = NACT_NU
C
      NSTNEO = K_NEO
C
C     --- MAY NEVER NEED 2-NUC DENSITY ---
C
      DOTWON =.TRUE.
      IF(NACT_NU.EQ.N1 .OR. (NA_NU+NB_NU.EQ.1)) THEN
         DOTWON = .FALSE.
      END IF
C
C     --- WORK OUT HOW MANY ELECTRONIC AND NUCLEAR DETERMINANTS
C                 AS WELL AS MIXED NUCLEAR-ELECTRONIC DETERMINANTS ---
C
      IF(SOME1) WRITE(IW,6090)
      CALL CNTDET(NALP,NBLP,NALP_NU,NBLP_NU,IDS,IIS,SOME1)
C
      NEDET = NCI
      NNDET = NCI_NU
C
      MIXDET = NCI*NCI_NU
C
      IF(SOME1) WRITE(IW,7000) MIXDET
C
      NEDET2 = (NEDET*NEDET+NEDET)/2
      NNDET2 = (NNDET*NNDET+NNDET)/2
      MIXDT2 = (MIXDET*MIXDET+MIXDET)/2
C
C     --- WORK OUT MAX SIZES FOR MEMORY ALLOCATION ---
C         (NA AND NA_NU ARE # OF ALPHA ACTIVE EL AND NU)
C
      MXNACT = MAX(NACT,NACT_NU)
      MXNTOT = MAX(NTOT,NTOT_NU)
      MXNTOT2 = MAX(NTOT2,NTOT2_NU)
      MXNA   = MAX(NA,NA_NU)
C
      NCI2 = (NCI*NCI+NCI)/2
      NCI2_NU = (NCI_NU*NCI_NU+NCI_NU)/2
      MXNCI2 = MAX(NCI2,NCI2_NU)
C
      M1 = NACT
      M2 = (M1*M1+M1)/2
      M4 = (M2*M2+M2)/2
C
      M1_NU = NACT_NU
      M2_NU = (M1_NU*M1_NU+M1_NU)/2
      M4_NU = (M2_NU*M2_NU+M2_NU)/2
C
      MXM1 = MAX(M1,M1_NU)
      MXM2 = MAX(M2,M2_NU)
      MXM4 = MAX(M4,M4_NU)
C
      M5 = M2*M2_NU
      MAX2IN = MAX(MXM4,M5)
      MXIA = MAX(MXM2,MXNCI2)
C
C     --- NOW ALLOCATE MEMORY ---
C                                                IDENTIFY POINTERS
      CALL VALFM(LOADFM)
      LDM1E = LOADFM + 1                    ! LDM1E: 1-EL DENSITY MATRIX
      LDM1N = LDM1E  + M2                   ! LDM1N: 1-NU DENSITY MATRIX
      LAST  = LDM1N  + M2_NU
      NEED = LAST - LOADFM - 1
C
      IECOUP    = LAST                      ! IECOUP: 1-EL COUPLING COEF.
      INCOUP    = IECOUP + NEDET2           ! INCOUP: 1-NU COUPLING COEF.
      LVEC      = INCOUP + NNDET2           ! LVEC  : NEO-CI EIGENVECTORS
      LEIG      = LVEC   + MIXDET*NSTNEO    ! LEIG  : NEO-CI EIGENVALUES
      LIA       = LEIG   + MIXDET
      IBO       = LIA    + MXIA/NWDVAR + 1  ! IBO   : MO SYMMETRY LABELS
      IACON1    = IBO    + MXNTOT           ! IACON1: ALPHA EL OCCS DET 1
      IBCON1    = IACON1 + MXNA             ! IBCON1: BETA EL OCCS  DET 1
      IACON2    = IBCON1 + MXNA             ! IACON2: ALPHA EL OCCS DET 2
      IBCON2    = IACON2 + MXNA             ! IBCON1: BETA EL OCCS  DET 2
      IACON1_NU = IBCON2 + MXNA             ! IACON1_NU: ALPHA NU OCCS DET 1
      IBCON1_NU = IACON1_NU + MXNA          ! IBCON1_NU: BETA NU OCCS DET 1
      IACON2_NU = IBCON1_NU + MXNA          ! IACON2_NU: ALPHA NU OCCS DET 2
      IBCON2_NU = IACON2_NU + MXNA          ! IBCON1_NU: BETA NU OCCS DET 2
      INDEX     = IBCON2_NU + MXNA          ! INDEX  : INDEXING ARRAY
      INDEXE    = INDEX  + ((MXNACT*(MXNACT+1))/2)**2
      INDEXN    = INDEXE + ((MXNACT*(MXNACT+1))/2)**2
      LAST      = INDEXN + ((MXNACT*(MXNACT+1))/2)**2
      NEED2 = LAST - LOADFM - NEED - 1

      LSINT1    = LAST                     ! LSINT1 : 1-PARTICLE MO INTS.
      LSINT2    = LSINT1 + MXM2            ! LSINT2 : 2-PARTICLE MO INTS.
      LXX       = LSINT2 + MAX2IN
      LIXX      = LXX    + NINTMX
      LHNMO     = LIXX   + NINTMX          ! LHNMO  : NUCLEAR EFF. FOCK MAT.
      LECIH     = LHNMO  + N2              ! LECIH  : ELEC. CI HAMILTONIAN
      LNCIH     = LECIH  + NEDET2          ! LNCIH  : NUCLEAR CI HAMILTONIAN
      LSCR      = LNCIH  + NNDET2          ! LSCR   : SCRATCH SPACE
      IWRK      = LSCR   + MIXDET*8        ! IWRK   : SCRATCH SPACE
      LAST      = IWRK   + MIXDET
      NEED3 = LAST - LOADFM - NEED - NEED2 - 1
C
C     --- ALLOCATE MEMORY FOR MIXED HAMILTONIAN SEPARATELY ---
C
      LNECIH    = LAST                     ! LNECIH : NEO-CI HAMILTONIAN
      LAST      = LNECIH + MIXDT2
      NEED4 = LAST - LOADFM - NEED - NEED2 - NEED3 - 1
C
C     --- CHECK IF IN MEMORY RUN POSSIBLE ---
C
      CALL GOTFM(NGOTMX)
      NEEDTL = NEED + NEED2 + NEED3 + NEED4
C
      IF(SOME1) THEN
         WRITE(IW,7025) NEEDTL
         WRITE(IW,7027) NGOTMX
      END IF
C
      IF(NEEDTL.GT.NGOTMX) THEN
         IF(MASWRK) WRITE(IW,7030)
         CALL ABRT
      ELSE
         IF(SOME1) WRITE(IW,7040)
      END IF
C
      CALL GETFM(NEED)
      CALL GETFM(NEED2)
      CALL GETFM(NEED3)
      CALL GETFM(NEED4)
C
C     --- FORM NUC-ELEC ORBITAL (NEO) CI HAMILTONIAN MATRIX ELEMENTS ---
C
      IF(NEOSCF.EQ.MCSCF) CLABEL = .TRUE.
      CALL DETFCI_RETMAT(CLABEL,X(LSINT1),X(LSINT2),X(LIA),
     *     X(LXX),X(LIXX),X(IACON1),X(IBCON1),X(IACON2),
     *     X(IBCON2),X(IACON1_NU),X(IBCON1_NU),X(IACON2_NU),
     *     X(IBCON2_NU),X(INDEX),X(LHNMO),X(LECIH),X(IECOUP),X(LNCIH),
     *     X(INCOUP),X(LNECIH),NALP,NBLP,NALP_NU,NBLP_NU,
     *     M1,M2,M4,M1_NU,M2_NU,M4_NU,M5,MXM1,MXM2,MAX2IN,NINTMX,
     *     MXNA,NCORSV,NACT,NA,NB,NCI,NACT_NU,
     *     NA_NU,NB_NU,NCI_NU)
C
      IF(SOME1) THEN
         WRITE(IW,7510)
         CALL TIMIT(1)
      END IF
C
C     --- ANY FROZEN CORE CONTRIBUTION PLUS CLASSICAL NUCLEAR ENERGY
C
      ECONST = ECORE + ENUCR
C
C     --- DIAGONALIZE MIXED ELECTRONIC-NUCLEAR CI MATRIX ---
C
      CALL GLDIAG(MIXDET,NSTNEO,MIXDET,X(LNECIH),X(LSCR),
     *            X(LEIG),X(LVEC),IGERR,X(IWRK))
C
      DO I=1,NSTNEO
         X(LEIG-1+I) = X(LEIG-1+I) + ECONST
         ESTATE(I) = X(LEIG-1+I)
      END DO
C
C     --- NO LONGER NEED CI HAMILTONIANS - RETURN MEMORY ---
C
      CALL RETFM(NEED4)
      CALL RETFM(NEED3)
C
C     --- PRINT OUT NEO-CI COEFFICENTS AND ASSOCIATED CONFIGS ---
C
      IF (NPASS.EQ.0 .OR. NPASS.EQ.-1 .AND. MASWRK) THEN
         IF(NPASS.EQ.0) THEN
            WRITE(IW,7073)
            WRITE(IW,7074)
         ELSE
            WRITE(IW,7083)
            WRITE(IW,7084)
         END IF
         CALL PRIDTN(X(LEIG),X(LVEC),NSTNEO,MIXDET,NACT,
     *               NA,NB,NACT_NU,NA_NU,NB_NU,PRTTOL_NEO,
     *               X(IACON1),X(IBCON1),X(IACON1_NU),X(IBCON1_NU))
         WRITE (IW,8000)
      END IF
C
C     --- WRITE NEO CI COEFFCIENTS TO DISK FILE ---
C
      CALL SEQOPN(NFT12,'CIVECTR','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQREW(NFT12)
      WRITE(NFT12) NSTNEO,MIXDET
      DO I=1,NSTNEO
         LCIVEC = LVEC + (I-1)*MIXDET
         CALL STFASE(X(LCIVEC),MIXDET,MIXDET,1)
         CALL SQWRIT(NFT12,X(LCIVEC),MIXDET)
      END DO
      CALL SEQREW(NFT12)
C
C     --- FORM NEO CI DENSITY MATRICES ---
C
C     1-PARTICLE DENSITIES PLUS OPTIONAL NATURAL ORBITALS
C                          AND/OR
C     1- AND 2-PARTICLE DENSITIES, PLUS MIXED EL-NU DENSITY,
C     PLUS STATE AVERAGING IF REQUESTED
C
      IF(NPASS.GT.0) GO TO 500
C
C     --- MEMORY FOR NATURAL ORBITALS ---
C
      CALL VALFM(LOADFM)
      LMOS  = LOADFM + 1
      LNOS  = LMOS   + MAX(L3,N3)
      LOCC  = LNOS   + MAX(L3,N3)
      IWRK  = LOCC   + MAX(L1,N1)
      LWRK  = IWRK   + MXM1      !MXM1 = MAX NUMBER OF ACTIVE ORBITALS
      LSCR  = LWRK   + MXM1*8
      LDAO  = LSCR   + MXM1
      IWRK2 = LDAO   + MAX(L2,N2)
      LSCR2 = IWRK2  + MAX(L1,N1)
      LAST  = LSCR2  + MAX(L1,N1)
      NEED6 = LAST - LOADFM - 1
      CALL GETFM(NEED6)
C
C     --- ONE-PARTICLE DENSITIES ---
C
      IRTBEG = 1
      IRTEND = K_NEO
C
      IF(SOME2) WRITE(IW,7550)
C
      DO 900 NROOT = IRTBEG, IRTEND
C
         IF(NFLGDM_NEO(NROOT) .EQ. 0) GO TO 900
C
C        --- FORM BOTH ELECTRONIC AND NUCLEAR 1-PARTICLE DENSITIES ---
C
         CALL VCLR(X(LDM1E),1,M2)
         CALL VCLR(X(LDM1N),1,M2_NU)
         NROOTTMP = NROOT
         CALL NEODM1(X(LVEC),X(LDM1E),X(LDM1N),
     *          X(INDEXE),X(INDEXN),X(IACON1),X(IBCON1),
     *          X(IACON2),X(IBCON2),NALP,NBLP,NACT,NA,NB,
     *          X(IACON1_NU),X(IBCON1_NU),X(IACON2_NU),X(IBCON2_NU),
     *          NALP_NU,NBLP_NU,NACT_NU,NA_NU,NB_NU,
     *          M1,M2,M1_NU,M2_NU,X(LIA),MXM2,MIXDET,
     *          NROOTTMP)
C
C        --- OUTPUT 1-PARTICLE DENSITIES ---
C
         IF(SOME2 .AND. NFLGDM_NEO(NROOT).EQ.2) THEN
            WRITE(IW,7600) NROOT
            CALL PRTRI(X(LDM1E),M1)
            WRITE(IW,7610) NROOT
            CALL PRTRI(X(LDM1N),M1_NU)
         END IF
C
C        --- ELECTRONIC NATURAL ORBITALS ---
C
         CALL DAREAD(IDAF,IODA,X(LMOS),L3,15,0)
         COROCC = TWO
         CALL NEONO(X(LDM1E),X(LMOS),X(LNOS),X(LOCC),X(IWRK),
     *              X(LWRK),X(LSCR),NCORSV,M1,M1,M2,L1,COROCC,IW)
C
C        --- WRITE 1-ELEC DENSTY AND NATURAL ORBITAL INFO FOR
C            IROOT_NEO TO DISK FOR USE ELSEWHERE ---
C
         IF(NROOT.EQ.IROOT_NEO .AND. NPASS.EQ.-1) THEN
            CALL DMTX(X(LDAO),X(LNOS),X(LOCC),NCORSV+M1,L1,L1)
            IF(MASWRK) THEN
               CALL DAWRIT(IDAF,IODA,X(LDAO),L2,16,0)
               CALL DAWRIT(IDAF,IODA,X(LNOS),L3,19,0)
               CALL DAWRIT(IDAF,IODA,X(LOCC),L1,21,0)
            END IF
         END IF
C
C        --- OUTPUT AND PUNCH TO .DAT ELECTRONIC NATURAL ORBITALS ---
C
         IF(SOME2) THEN
            WRITE(IW,7620) NROOT
            CALL PREV(X(LNOS),X(LOCC),L0,L1,L1)
         END IF
         IF(SOME2) THEN
            CALL TMDATE(TIMSTR)
            WRITE(IP,7630) TIMSTR,TITLE,NROOT,X(LEIG+NROOT-1)
            WRITE(IP, FMT='(8H $VEC   )')
            CALL PUSQL(X(LNOS),L0,L1,L1)
            WRITE(IP, FMT='(5H $END)')
         END IF
C
C        --- NUCLEAR NATURAL ORBITALS ---
C
         CALL DAREAD(IDAF,IODA,X(LMOS),N3,444,0)
         COROCC = ONE
         CALL NEONO(X(LDM1N),X(LMOS),X(LNOS),X(LOCC),X(IWRK),
     *              X(LWRK),X(LSCR),NCORSV_NU,M1_NU,M1_NU,
     *              M2_NU,N1,COROCC,IW)
C
C        --- OUTPUT AND PUNCH TO .DAT NUCLEAR NATURAL ORBITALS ---
C
         IF(SOME2) THEN
            WRITE(IW,7650) NROOT
            CALL PREVNU(X(LNOS),X(LOCC),N1,N1,L1+1,L1+N1)
         END IF
         IF(SOME2) THEN
            CALL TMDATE(TIMSTR)
            WRITE(IP,7660) TIMSTR,TITLE,NROOT,X(LEIG+NROOT-1)
            WRITE(IP, FMT='(8H $VECPRO)')
            CALL PUSQL(X(LNOS),N1,N1,N1)
            WRITE(IP, FMT='(5H $END)')
         END IF
C
  900 CONTINUE
C
C     --- DONE WITH 1-PARTICLE DENSITIES AND NATURAL ORBITALS ---
C
      IF(SOME2) WRITE(IW,7950)
C
      IF(NPASS.EQ.-1 .AND. NEOSCF.EQ.MCSCF) THEN
C
C        --- OUTPUT NEO-MCSCF OPTIMIZED ORBITALS HERE ---
C
C        -REUSE SOME MEMORY
         LS = LDAO
         LQ = LNOS
         LEIG = LOCC
C
C        --- NO CANONICALIZATION DONE AS OF YET ---
C
C        -NO NEO-MCSCF SYMMETRY YET, EVEN SO GO THROUGH THE MOTIONS
         CALL DAREAD(IDAF,IODA,X(LS),L2,12,0)
         CALL DAREAD(IDAF,IODA,X(LMOS),L3,15,0)
         CALL DAREAD(IDAF,IODA,X(LQ),L3,45,0)
         CALL SYMMOS(X(IWRK2),X(LQ),X(LS),X(LMOS),X(LSCR2),L0,L1,L0,L1)
         CALL DAWRIT(IDAF,IODA,X(IWRK2),L1,255,1)
C
         CALL VCLR(X(LEIG),1,L1)
         IF(MASWRK) THEN
            WRITE(IW,7960)
            WRITE(IW,7965)
            CALL PREVS(X(LMOS),X(LEIG),X(IWRK2),L0,L1,L1,1)
         END IF
         IF(MASWRK) THEN
            CALL TMDATE(TIMSTR)
            WRITE(IP,7970) TIMSTR
            WRITE(IP,7975) TITLE,ETOT,ENUCR
            WRITE(IP,FMT='(8H $VEC   )')
            CALL PUSQL(X(LMOS),L0,L1,L1)
            WRITE(IP, FMT='(5H $END)')
         END IF
C
C        --- NOTE: GIVEN THE CURRENT RESTRICTIONS REGARDING ALL NUC
C            MOS BEING ACTIVE, NUCLEAR MOS ARE THE SAME AS THOSE READ IN.
C            STILL OUPUT USEFUL FOR USERS RECORD ---
C
         CALL DAREAD(IDAF,IODA,X(LMOS),N3,444,0)
         CALL VCLR(X(LEIG),1,N1)
         IF(MASWRK) THEN
            WRITE(IW,7980)
            WRITE(IW,7965)
            WRITE(IW,7985)
            CALL PREVNU(X(LMOS),X(LEIG),N1,N1,L1+1,L1+N1)
         END IF
         IF(MASWRK) THEN
            CALL TMDATE(TIMSTR)
            WRITE(IP,7970) TIMSTR
            WRITE(IP,7975) TITLE,ETOT,ENUCR
            WRITE(IP,FMT='(8H $VECPRO)')
            CALL PUSQL(X(LMOS),N1,N1,N1)
            WRITE(IP, FMT='(5H $END)')
         END IF
C
      END IF
C
C     --- RETURN SOME MEMORY ---
C
      CALL RETFM(NEED6)
C
  500 CONTINUE
C
C     --- PUT ENERGY IN COMMON BLOCK FOR POSSIBLE MCSCF OR
C         NUMERICAL GRADIENT. STATE-AVERGING OF THE ENERGY IS
C         CARRIED OUT HERE IF REQUIRED ---
C
C     --- NO MCSCF, PUT CHOSEN ROOT'S ENERGY IN COMMON BLOCKS --
C
      ETOT = ESTATE(IROOT_NEO)
      E = ETOT
      EELCT = ETOT - ENUCR
C
C     --- MCSCF, DO POSSIBLE STATE-AVERAGING OF THE ENERGY ---
C
      IF(NPASS.NE.-1 .AND. NEOSCF.EQ.MCSCF) THEN
C
         MXSTAT=0
         MXNZW=0
         DO 200 I=MXRT,1,-1
            IF(WSTATE_NEO(I).GT.ZERO) THEN
               IF(MXSTAT.EQ.0) MXSTAT=I
               MXNZW=MXNZW+1
            END IF
  200    CONTINUE
         IF(SOME1) WRITE(IW,8520) MXNZW
C
         E = ZERO
         NXTR = 0
         DO 300 IST=1,K_NEO
            NXTR = NXTR + 1
            IF(WSTATE_NEO(NXTR).GT.ZERO) THEN
               E = E + WSTATE_NEO(NXTR) * ESTATE(IST)
               IF(SOME1) WRITE(IW,8525) IST,ESTATE(IST),
     *                                 WSTATE_NEO(NXTR)
            END IF
  300    CONTINUE
C
         ETOT = E
         EELCT = ETOT - ENUCR
C
         IF(NXTR.LT.MXSTAT) THEN
            IF(MASWRK) WRITE(IW,8530) NXTR,MXSTAT
            CALL ABRT
         END IF
C
      END IF
C
C     --- FOR NEO-MCSCF RUNS WILL NEED 2-PARTICLE DENSITIES ---
C
      IF(NEOSCF.EQ.MCSCF .AND. NPASS.NE.-1) THEN
C
C        --- ADDITIONAL MEMORY NEEDED ---
C
         CALL VALFM(LOADFM)
         LDMX   = LOADFM + 1
         LEDM1  = LDMX   + M2*M2_NU
         LEDM2  = LEDM1  + M2
         LNDM1  = LEDM2  + M4
         LNDM2  = LNDM1  + M2_NU
         LSDMX  = LNDM2  + M4_NU
         LSEDM1 = LSDMX  + M2*M2_NU
         LSEDM2 = LSEDM1 + M2
         LSNDM1 = LSEDM2 + M4
         LSNDM2 = LSNDM1 + M2_NU
         LXX    = LSNDM2 + M4_NU
         LIXX   = LXX    + NINTMX
         LDNWRK = LIXX   + NINTMX
         LLABMO = LDNWRK + MXNTOT2
         LLBABL = LLABMO + MAX(L1,N1)
         LLBIRP = LLBABL + MAX(M1,M1_NU)
         LSYIRP = LLBIRP + 12
         LAST   = LSYIRP + 12
         NEED7  = LAST - LOADFM - 1
         CALL GETFM(NEED7)
C
C        --- 1- AND 2-PARTICLE DENSITIES AND STATE-AVERAGING  ---
C
         CALL NEODM2(X(IECOUP),X(INCOUP),X(LVEC),X(LDMX),X(LEDM1),
     *        X(LEDM2),X(LNDM1),X(LNDM2),X(INDEX),X(INDEXE),X(INDEXN),
     *        X(IACON1),X(IBCON1),X(IACON2),X(IBCON2),
     *        NALP,NBLP,NACT,NA,NB,
     *        X(IACON1_NU),X(IBCON1_NU),X(IACON2_NU),X(IBCON2_NU),
     *        NALP_NU,NBLP_NU,NACT_NU,NA_NU,NB_NU,M1,M2,M4,
     *        M1_NU,M2_NU,M4_NU,X(LIA),MXM2,MIXDET,K_NEO,
     *        WSTATE_NEO,X(LSDMX),X(LSEDM1),X(LSEDM2),X(LSNDM1),
     *        X(LSNDM2),DOTWON)
C
C        --- WRITE ELECTRONIC DENSITY MATRICES TO DISK ---
C
         CALL DETGRP(GRPDET,X(LLABMO),X(LLBABL),PTGRP,X(LLBIRP),
     *               X(LSYIRP),NSYM,NIRRP,L1,NACT,NCORSV)
C
CSPW     CALL DAWRIT(IDAF,IODA,X(LSEDM1),M2,454,0)
CSPW     CALL DAWRIT(IDAF,IODA,X(LSEDM2),M4,455,0)
C
         CUTOFF = MAX(1.0D-11,10.0D+00**(-ICUT))
         IF(SOME1) WRITE(IW,9370) X(LSYIRP),PTGRP
C
         CALL SEQOPN(NFT15,'WORK15','UNKNOWN',.FALSE.,'UNFORMATTED')
         CALL SEQREW(NFT15)
         CALL WTDNEO(EXETYP,X(LSEDM1),X(LSEDM2),X(LLBABL),X(LXX),
     *               X(LIXX),NINTMX,LABSIZ,M1,M2,M4,
     *               X(LDNWRK),NTOT2,NCORSV,CUTOFF,NFT15,NRECO,NDM2O)
         CALL SEQREW(NFT15)
C
         IF(SOME1) WRITE(IW,9380) NDM2O,NRECO,NFT15
         IF(SOME1) WRITE(IW,9390)
C
C        --- WRITE NUCLEAR DENSITY MATRICES TO DISK ---
C
CSPW    2-NUCLEUS DENSITY SHOULD PROBABLY HAVE REGULAR DISK FILE (NOT DAF)
         CALL DAWRIT(IDAF,IODA,X(LSNDM1),M2_NU,456,0)
         IF(DOTWON) CALL DAWRIT(IDAF,IODA,X(LSNDM2),M4_NU,457,0)
C
C        --- WRITE ELECTRONIC-NUCLEAR DENSITY MATRICES ---
C
CSPW    MIXED DENSITY SHOULD PROBABLY HAVE REGULAR DISK FILE (NOT DAF)
         CALL DAWRIT(IDAF,IODA,X(LSDMX),M5,458,0)
C
         CALL RETFM(NEED7)
C
      END IF
C
C     --- RETURN REST OF MEMORY ---
C
      CALL RETFM(NEED2)
      CALL RETFM(NEED)
C
      RETURN
 6090 FORMAT(/,5X,45('-'),
     *       /,5X,'NUMBER OF ELECTRONIC AND NUCLEAR DETERMINANTS',
     *       /,5X,45('-'))
 7000 FORMAT(/,1X,62('-'),/1X,'TOTAL NUMBER OF NUCLEAR-ELECTRONIC ',
     *       ' DETERMINANTS IS ',I10)
 7025 FORMAT(/,1X,'IN MEMORY NEO DETERMINANT FULL CI REQUIRES',I14,
     *       ' WORDS')
 7027 FORMAT(1X,'TOTAL MEMORY AVAILABLE IS',17X,I14,' WORDS',
     *       /,1X,62('-'))
 7030 FORMAT(/,1X,'NEO-FCI REQUIRES MORE MEMORY - ABORTING')
 7040 FORMAT(/,1X,'NEO-FCI HAMILTONIAN WILL BE FORMED IN MEMORY ...',/)
 7073 FORMAT(/,1X,'NEO-FCI EIGENVECTORS AND EIGENVALUES')
 7074 FORMAT(1X,36('-'),/)
 7083 FORMAT(/,1X,'NEO-FCI BASED ON OPTIMIZED ORBITALS')
 7084 FORMAT(1X,35('-'),/)
 7510 FORMAT(1X,'DONE WITH NEO DETERMINANT CI HAMILTONIAN ',
     *       'COMPUTATION .....'/)
 7550 FORMAT(5X,31(1H-)/
     * 5X,'FORM NEO ONE-PARTICLE DENSITIES'/
     * 5X,31(1H-))
 7600 FORMAT(/,1X,'ONE-ELECTRON DENSITY (MO BASIS) FOR NEO STATE',I4)
 7610 FORMAT(/,1X,'ONE-NUCLEUS DENSITY (MO BASIS) FOR NEO STATE',I4)
 7620 FORMAT(/,1X,'ELECTRONIC NATURAL ORBITALS FOR NEO STATE',I4)
 7630 FORMAT('--- ELECTRONIC NEO-CI NATURAL ORBITALS --- GENERATED AT '
     *       ,3A8/10A8/'STATE ',I2,' E=',F20.10)
 7650 FORMAT(/,1X,'NUCLEAR NATURAL ORBITALS FOR NEO STATE',I4)
 7660 FORMAT('--- NUC NEO-CI NATURAL ORBITALS --- GENERATED AT ',
     *       3A8/10A8/'STATE ',I2,' E=',F20.10)
 7950 FORMAT(/1X,'..... DONE WITH NEO ONE-PARTICLE DENSITIES ..... '/)
 7960 FORMAT(/10X,39(1H-)/10X,'NEO-MCSCF OPTIMIZED ELECTRONIC ORBITALS'/
     *        10X,39(1H-) )
 7965 FORMAT(/10X,'THESE ORBITALS ARE NOT CANONICALIZED ...')
 7970 FORMAT('--- OPTIMIZED MCSCF MO-S --- GENERATED AT ',3A8)
 7975 FORMAT(10A8/'E(MCSCF)=',F20.10,', E(NUC)=',F20.10)
 7980 FORMAT(/10X,36(1H-)/10X,'NEO-MCSCF OPTIMIZED NUCLEAR ORBITALS'/
     *        10X,36(1H-) )
 7985 FORMAT(/10X,'AS ALL NMOS ARE ACTIVE THEY ARE THE SAME AS THOSE',
     *            ' READ IN ...')
C
 8000 FORMAT(/,1X,'.... DONE WITH NEO-FCI ....',/)
 8520 FORMAT(1X,'STATE AVERAGING WILL OCCUR OVER',I4,' ROOT(S)')
 8525 FORMAT(1X,'STATE=',I4,'   ENERGY=',F20.10,'   WEIGHT=',F8.5)
 8530 FORMAT(/1X,'***** ERROR *****'/
     *       1X,'THIS RUN FOUND',I5,' CI EIGENVECTORS,'/
     *       1X,'BUT YOU REQUESTED STATE AVERAGING OF',I5,' ROOTS.'/
     *       1X,'PLEASE EXAMINE YOUR CHOICE OF -NSTATE- IN $NEODET.'/)
 9370 FORMAT(/1X,'SIEVING THE ',A4,
     *          ' SYMMETRY NONZERO DENSITY ELEMENTS IN GROUP ',A8)
 9380 FORMAT(1X,I10,' NONZERO DM2 ELEMENTS WRITTEN IN',I8,
     *          ' RECORDS TO FILE',I3)
 9390 FORMAT(1X,'..... DONE WITH 1 AND 2 PARTICLE DENSITY MATRIX .....')
      END
C
C*MODULE NEOFCI  *DECK NDETIN
      SUBROUTINE NDETIN(NPRINT,GPNAME)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DOUBLE PRECISION NEOSCF,MCSCF,NEOCI
C
      LOGICAL SOME,PURES,GOPARR,DSKWRK,MASWRK,ABEL,DAVID
C
      PARAMETER (MXATM=2000, MXRT=100, MXAO=8192, MXSH=5000)
C
      COMMON /DAVSON/ DAVID
      COMMON /DETNEO/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,S,SZ,
     *                GRPDET,STSYM,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORB,
     *                NA,NB,K,KST,IROOT,IPURES,MAXW1,NITER,MAXP,NCI,
     *                IGPDET,KSTSYM
      COMMON /FMCOM / X(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,MA,MB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IJKT,IDAF,NAV,IODA(950)
      COMMON /NEOCIR/ IEROOT
      COMMON /NEOMTD/ NEOSCF,NEOCI
      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMMOL/ GROUP,COMPLEX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
      DATA MCSCF/8HMCSCF   /
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
C
      PARAMETER (NNAM=19)
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DIMENSION FANT(8),LFANT(8),GANT(27),LGANT(8)
C
C     --- READ $NEODET USER INPUT GROUP - THIS INCLUDES NUMBER OF
C         ACTIVE NUCLEI, NUCLEAR ORBITAL ACTIVE ORBITALS, AND THE
C         NUCLEAR SPIN, AS WELL AS NUMBER OF NEO STATES,
C         CONTROL OF DAVIDSON DIAG. (FUTURE), CONTROL OF
C         NEO STATE AVERAGING, AND CONTROL OF CALCULATION OF NEO
C         NATURAL ORBITALS.
C
C         $NEODET REQUESTS FOR STATE-AVERAGING, NUMBER OF STATES,
C         AND ONE-PARTICLE PROPERTIES WILL OVERIDE ANY $DET
C         RELATED INPUT ---
C
C         THIS ROUTINE IS DETINP WITH VERY MINOR MODIFICATIONS - SPW
C
      DATA QNAM/8HGROUP   ,8HISTSYM  ,
     *          8HNCORE   ,8HNACT    ,8HNUCLEI  ,8HSZ      ,
     *          8HNSTATE  ,8HNSTGSS  ,8HNHGSS   ,8HMXXPAN  ,
     *          8HITERMX  ,8HCVGTOL  ,8HPRTTOL  ,
     *          8HIROOT   ,8HNFLGDM  ,8HPURES   ,8HWSTATE  ,
     *          8HDAVID   ,8HIEROOT  /
      DATA KQNAM/5,1,1,1,1,3,  1,1,1,1,1,3,3,  1,-1,0,-3,0,1/
C
      DATA HESS/8HHESSIAN /
      DATA FANT/8HC1      ,8HCI      ,8HCS      ,8HC2      ,
     *          8HD2      ,8HC2V     ,8HC2H     ,8HD2H     /
      DATA LFANT/1,1,1,1,2,2,2,3/
      DATA GANT/8HA       ,8HAG      ,8HAU      ,8HA'      ,
     *          8HA"      ,8HA       ,8HB       ,8HA       ,
     *          8HB1      ,8HB2      ,8HB3      ,8HA1      ,
     *          8HA2      ,8HB1      ,8HB2      ,8HAG      ,
     *          8HBG      ,8HBU      ,8HAU      ,8HAG      ,
     *          8HB1G     ,8HB2G     ,8HB3G     ,8HAU      ,
     *          8HB1U     ,8HB2U     ,8HB3U     /
      DATA LGANT/0,1,3,5,7,11,15,19/
C
      SOME = MASWRK  .AND.  NPRINT.NE.-5  .AND.  NPRINT.NE.-23
C
      IF(SOME) WRITE(IW,9000)
C
C          SET UP INPUT TO SPECIFY THE CI SPACE
C          SELECT POINT GROUP AND NUCLEI/ORBITAL COUNTS
C
      GRPDET = FANT(1)
      IF(IGROUP.EQ.1)                GRPDET = FANT(1)
      IF(IGROUP.EQ.3)                GRPDET = FANT(2)
      IF(IGROUP.EQ.2)                GRPDET = FANT(3)
      IF(IGROUP.EQ.4.AND.NAXIS.EQ.2) GRPDET = FANT(4)
      IF(IGROUP.EQ.8.AND.NAXIS.EQ.2) GRPDET = FANT(5)
      IF(IGROUP.EQ.7.AND.NAXIS.EQ.2) GRPDET = FANT(6)
      IF(IGROUP.EQ.6.AND.NAXIS.EQ.2) GRPDET = FANT(7)
      IF(IGROUP.EQ.9.AND.NAXIS.EQ.2) GRPDET = FANT(8)
      IF (RUNTYP.EQ.HESS  .OR.  NT.EQ.1) GRPDET=FANT(1)
      KSTSYM = 1
C
      NCORE  = 0
      NACT   = 0
      NUCLEI = 0
      SZ     = (NUMULT-1)/TWO
C
C          SET UP INPUT TO CONTROL THE DIAGONALIZATION
C
      NSTATE = 1
      NSTGSS = 1
      NHGSS  = 300
      MXXPAN = 10
      ITERMX = 100
      CVGTOL = 1.0D-05
      PRTTOL = 0.05D+00
      DAVID  = .FALSE.
C
C          SET UP INPUT TO CONTROL THE FIRST ORDER DENSITY COMPUTATION
C
      IROOT=1
      IEROOT=0
      KQNAM(15)=MXRT*10 + 1
      DO 5 I=1,MXRT
         NFLGDM(I) = 0
    5 CONTINUE
      NFLGDM(1)=1
C
C          SET UP INPUT TO CONTROL THE SECOND ORDER DENSITY COMPUTATION
C
      PURES = .TRUE.
      KQNAM(17)=MXRT*10 + 3
      CALL VCLR(WSTATE,1,MXRT)
      WSTATE(1) = ONE
C
      CALL NAMEIO(IR,JRET,GPNAME,NNAM,QNAM,KQNAM,
     *            GRPDET,KSTSYM,NCORE,NACT,NUCLEI,SZ,NSTATE,NSTGSS,
     *            NHGSS,MXXPAN,ITERMX,CVGTOL,PRTTOL,IROOT,NFLGDM,
     *            PURES,WSTATE,DAVID,IEROOT,
     *    0,0,0,0,0,
     *    0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,
     *    0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0)
C
      IF(JRET.EQ.2) THEN
         IF(MASWRK) WRITE(IW,9010) GPNAME
         CALL ABRT
      END IF
C
C     --- NEO CANNOT USE SYMMETRY YET - OVERRIDE ANY THE USER HAS INPUT ---
      GRPDET = FANT(1)
      KSTSYM = 1
C
      IGPDET = -1
C
C        THE INPUT FOR C2H IS SUPPOSED TO BE IDENTICAL TO THE GUGA
C        ORDER, NAMELY 1,2,3,4=AG,BU,BG,AU, BUT THE CI CODE WANTS
C        THE ORDER OF  1,2,3,4=AG,BG,BU,AU.  SEE ALSO GAJASW ROUTINE.
C
      IF (GRPDET.EQ.FANT(7)) THEN
         MODI = KSTSYM
         IF(KSTSYM.EQ.2) MODI=3
         IF(KSTSYM.EQ.3) MODI=2
         KSTSYM=MODI
      END IF
      DO I=1,8
         IF (GRPDET.EQ.FANT(I)) THEN
            IGPDET=LFANT(I)
            STSYM = GANT(LGANT(I)+KSTSYM)
         ENDIF
      END DO
      IF (IGPDET.EQ.-1) THEN
         IF(MASWRK) WRITE(IW,*) '$NEODET POINT GROUP IS UNRECOGNIZED!'
         CALL ABRT
      ENDIF
      IF (GRPDET.EQ.FANT(1).AND.KSTSYM.GT.1) THEN
         IF(MASWRK) WRITE(IW,*) '$NEODET STATE SYMM IS NOT CORRECT',
     *                          ' IRREP'
         CALL ABRT
      ENDIF
      IF (KSTSYM.GT.(2**IGPDET)) THEN
         IF(MASWRK) WRITE(IW,*)
     *       '$NEODET STATE SYMMETRY IS TOO LARGE FOR THIS GROUP'
         CALL ABRT
      ENDIF
C
C READ MO SYMMETRIES AND WRITE TO DIRECT ACCESS FILE,
C BUT DON'T DO IF THIS IF WE ARE PEEKING AT THE $DET INPUT
C AT THE MOST EARLY STAGES OF AN INITIAL MCSCF RUN.
C
      IF (NPRINT.EQ.-23) GOTO 1314
C
CSPW  N0 = NQMT
      N0 = NUMNB
      N1 = NUMNB
      N2 = (N1*N1+N1)/2
      N3 =  N1*N1
C
      CALL VALFM(LOADFM)
      LMOLAB = LOADFM + 1
      LMOIRP = LMOLAB + N1
      LVEC   = LMOIRP + N1
      LS     = LVEC   + N3
      LQ     = LS     + N2
      LWRK   = LQ     + N3
      LMODEG = LWRK   + N1
      LAST   = LMODEG + N1
      NEEDD   = LAST - LOADFM - 1
      CALL GETFM(NEEDD)
C
      IF(GRPDET.EQ.FANT(1)) THEN
         CALL C1DET(X(LMOIRP),X(LMOLAB),N0)
      ELSE
         CALL DAREAD(IDAF,IODA,X(LVEC),N3,444,0)
         CALL DAREAD(IDAF,IODA,X(LS),N2,441,0)
         CALL DAREAD(IDAF,IODA,X(LQ),N3,443,0)
         CALL TRFSYM(X(LMOLAB),X(LMOIRP),X(LMODEG),X(LQ),X(LS),X(LVEC),
     *               X(LWRK),IA,N0,N1,N0,N1)
      END IF
C
C  CHANGE ORBITAL SYMMETRY LABELS FROM GAMESS TO JAKAL VALUES
CSPW NEED TO SET UP FOR NEO WHEN SYMMETRY IS INCLUDED
C
C     DO I=1,NUMNB
C        X(LMOIRP-1+I) = 1
C     END DO
C     DO I=1,N1
C     WRITE(6,*)'LMOIRP=',X(LMOIRP-1+I)
C     END DO
      CALL PUTINT(X(LMOIRP),N1)
C     CALL GAJASW(X(LMOIRP),NUMNB,GRPDET)
C
      CALL DAWRIT(IDAF,IODA,X(LMOIRP),N1,453,1)
C
C     1.  SET NUCCSV,NCOR,NUORBA,NORB,NA,NB FOR DETERMINANT SPECIFICATION
C         CHECK INPUT, AND COPY INTO INTERNALLY USED VARIABLE NAMES.
C         NCOR WILL BE SET TO ZERO TO DROP CORES, SO NUCCSV SAVES THIS.
C
 1314 CONTINUE
C
      IF(NPRINT.NE.-23  .AND.
     *       (NCORE.LT.0  .OR.  NACT.LE.0  .OR.  NUCLEI.LE.0)) THEN
         IF(MASWRK) WRITE(IW,9020) GPNAME,NCORE,NACT,NUCLEI
         CALL ABRT
      END IF
      NCORSV = NCORE
      NCOR   = NCORE
      NORB   = NCORE + NACT
      NHIGH = INT(SZ+SZ+0.0001D+00)
      NB = (NUCLEI-NHIGH)/2
      NA = NB+NHIGH
      MA = NA+NCORSV
      MB = NB+NCORSV
      NUCTOT = 2*NCOR+NA+NB
      NERR=0
      IF(NUCTOT.NE.NQMNUC)  NERR=1
      IF(NUCLEI.NE.NA+NB) NERR=1
      IF(NA.LT.NB)      NERR=1
      IF(NA.LE.0)       NERR=1
      IF(NB.LT.0)       NERR=1
      IF(NPRINT.NE.-23  .AND.  NERR.GT.0) THEN
         IF(MASWRK) WRITE(IW,9030) NCORE,NUCLEI,SZ,NUMULT
         CALL ABRT
      END IF
      S = (NUMULT-1)/TWO
C
C        2. SET K,KST,MAXW1,NITER,MAXP,CRIT,PRTTOL FOR DIAGONALIZATION
C
      K     = NSTATE
      KST   = MAX(NSTGSS,K)
      MAXP  = MAX(MXXPAN,2*KST)
      MAXW1 = NHGSS
      NITER = ITERMX
      CRIT  = CVGTOL
C
      IF(NPRINT.NE.-23  .AND.  K.GT.MXRT) THEN
         IF(MASWRK) WRITE(IW,9035) K,MXRT
         CALL ABRT
         STOP
      END IF
C
C        3. SETUP FOR 1ST ORDER DENSITY COMPUTATION
C
      IF(IEROOT.GT.0. .AND. NEOSCF.EQ.MCSCF) THEN
         IROOT = IEROOT
      ELSE
         IEROOT = 0
      END IF
      IF(IROOT.GT.NSTATE) THEN
         IF(MASWRK) WRITE(IW,9036) IROOT,NSTATE
         CALL ABRT
      END IF
      IF(IROOT.GT.MXRT) THEN
         IF(MASWRK) WRITE(IW,9037) IROOT,MXRT
         CALL ABRT
      END IF
      IF(NFLGDM(IROOT).EQ.0) NFLGDM(IROOT)=1
C
C        4. SETUP FOR STATE-AVERAGING 1ST AND 2ND ORDER DENSITIES.
C        -IWTS- INDEXES THE NON-ZERO ELEMENTS OF -WSTNUC-
C
      IPURES=0
      IF(PURES) IPURES=1
      MXSTAT=0
      WSUM = ZERO
      DO 15 I=1,MXRT
         IF(WSTATE(I).GT.ZERO) THEN
            IF(I.LE.NSTATE) THEN
               WSUM = WSUM + WSTATE(I)
               MXSTAT = MXSTAT+1
               IWTS(MXSTAT) = I
            ELSE
               IF(MASWRK) WRITE(IW,9040) NSTATE
               CALL ABRT
            END IF
         END IF
         IF(WSTATE(I).LT.ZERO) THEN
            IF(MASWRK) WRITE(IW,9050)
            CALL ABRT
         END IF
   15 CONTINUE
      SCALE = ONE/WSUM
      CALL DSCAL(MXRT,SCALE,WSTATE,1)
C
C        IF RUNNING SILENTLY, RETURN WITHOUT PRINTING ANYTHING
C
      IF(NPRINT.EQ.-23) RETURN
C
      IF(NPROC.GT.1) THEN
         IF(MASWRK) WRITE(IW,9070)
         CALL ABRT
         STOP
      END IF
C
C     THE COMPUTATION OF THE ENERGY GRADIENT REQUIRES A SYMMETRIC
C     LAGRANGIAN (SINCE GAMESS DOES NOT AT PRESENT DO THE CPHF
C     PROBLEM TO OBTAIN ORBITAL PERTURBATIONS).  ALTHOUGH A STATE
C     AVERAGED MCSCF RUN RESULTS IN A SYMMETRIC "AVERAGE LAGRANGIAN",
C     THE LAGRANGIAN FOR ANY SINGLE STATE USING THE AVERAGED ORBITALS
C     WILL *NOT* BE SYMMETRIC.  HENCE ANY JOB THAT ATTEMPTS TO BOTH
C     STATE AVERAGE AND DO A GRADIENT SHOULD BE FLUSHED.
C
      CALL DERCHK(MAXDER)
      IF(MXSTAT.GT.1  .AND.  MAXDER.GT.0) THEN
         IF (SOME) WRITE(IW,9060) RUNTYP,MXSTAT
         CALL ABRT
         STOP
      END IF
C
      IF(SOME) THEN
         WRITE(IW,9100) GRPDET,STSYM,NCOR,NACT,
     *                  NA+NCOR,NA,NB+NCOR,NB,NORB
C
         WRITE(IW,9160) NCOR,NACT
         CALL MOSYPR(X(LMOLAB),NCOR,NACT)
C
         WRITE(IW,9108)
CSPW     WHEN DAVIDSON IMPLEMENTED THE FOLLOWING LINE WILL BE NEEDED
CSPW     WRITE(IW,9110) K,KST,MAXP,MAXW1,NITER,CRIT
         WRITE(IW,9110) K
         WRITE(IW,9120) IROOT
         WRITE(IW,9130) (NFLGDM(II),II=1,K)
         IF(NEOSCF.EQ.MCSCF .AND. IEROOT.GT.0) THEN
CSPW        WRITE(IW,9140) PURES
            WRITE(IW,9150) (IWTS(II),WSTATE(IWTS(II)),II=1,MXSTAT)
            IF(MXSTAT.GT.1) WRITE(IW,9155) IEROOT
         END IF
      END IF
C
C     --- CURRENTLY NEO-MCSCF REQUIRES NUC ACTIVE SPACE SPANS ALL MOS ---
C
      IF(NEOSCF.EQ.MCSCF .AND. NACT.NE.NUMNB) THEN
         WRITE(IW,9200)
         CALL ABRT
      END IF
C
      CALL RETFM(NEEDD)
      RETURN
C
 9000 FORMAT(/5X,20(1H-)/
     *       5X,'NUCLEAR ACTIVE SPACE'/
     *       5X,20(1H-))
 9010 FORMAT(/1X,'**** ERROR, THIS RUN REQUIRES INPUT OF A $',A8,
     *          'GROUP')
 9020 FORMAT(/1X,'**** ERROR, THIS RUN DOES NOT CORRECTLY SPECIFY',
     *          ' THE FULL NUCLEAR CI SPACE'/
     * 1X,'CHECK $',A8,' INPUT: NCORE=',I4,' NACT=',I4,' NUCLEI=',I4)
 9030 FORMAT(/1X,'**** ERROR, $NUCDET INPUT NCORE=',I4,' NUCLEI=',I4,
     *          ' SZ=',F6.3/
     *       1X,' IS INCONSISTENT WITH $NEO INPUT NUMULT=',I4)
 9035 FORMAT(/1X,'***** ERROR, REQUESTED NUMBER OF CI ROOTS=',I5/
     *        1X,'EXCEEDS THE DIMENSION LIMIT FOR NUMBER OF STATES',I5)
 9036 FORMAT(/1X,'**** ERROR, YOUR STATE SELECTED FOR PROPERTIES=',I5/
     *        1X,'EXCEEDS THE NUMBER OF ROOTS YOU REQUESTED=',I5)
 9037 FORMAT(/1X,'**** ERROR, YOUR STATE SELECTED FOR PROPERTIES=',I5/
     *        1X,'EXCEEDS THE DIMENSION LIMIT FOR NUMBER OF STATES',I5)
 9040 FORMAT(/1X,'**** ERROR, WEIGHTS ASSIGNED TO STATES HIGHER',
     *          ' THAN NSTNEO=',I5)
 9050 FORMAT(/1X,'**** ERROR, NEGATIVE VALUE FOR -WSTNUC- ???')
 9060 FORMAT(/1X,'**** ERROR, RUNTYP=',A8,' REQUIRES ENERGY GRADIENT.'/
     *       1X,'THIS IS IMPOSSIBLE WHILE STATE AVERAGING. NAVG=',I5)
 9070 FORMAT(/1X,'**** ERROR, DETERMINANT CODE DOES NOT IN PARALLEL.'/
     *       1X,'PLEASE RUN USING $MCSCF CISTEP=GUGA, AND A $DRT.')
C
 9100 FORMAT(/1X,'THE POINT GROUP                     =',3X,A8/
     *       1X,'THE STATE SYMMETRY                  =',3X,A8/
     *       1X,'NUMBER OF CORE NUCLEAR ORBITALS     =',I5/
     *       1X,'NUMBER OF ACTIVE NUCLEAR ORBITALS   =',I5/
     *       1X,'NUMBER OF ALPHA NUCLEI              =',I5,
     *          ' (',I4,' ACTIVE)'/
     *       1X,'NUMBER OF BETA NUCLEI               =',I5,
     *          ' (',I4,' ACTIVE)'/
     *       1X,'NUMBER OF OCCUPIED NUCLEAR ORBITALS =',I5)
 9108 FORMAT(/5X,38(1H-)/
     *       5X,'NEO-FCI DIAGONALIZATION AND PROPERTIES'/
     *       5X,38(1H-)/)
 9110 FORMAT(1X,'NUMBER OF NEO-CI STATES REQUESTED       =',I5)
C9110 FORMAT(1X,'NUMBER OF NEO-CI STATES REQUESTED       =',I5/
C    *       1X,'NUMBER OF NEO-CI STARTING VECTORS       =',I5/
C    *       1X,'MAX. NO. OF NEO-CI EXPANSION VECTORS    =',I5/
C    *       1X,'SIZE OF INITIAL NEO-CI GUESS MATRIX     =',I5/
C    *       1X,'MAX. NO. OF NEO-CI ITERS/STATE          =',I5/
C    *       1X,'NEO-CI DIAGONALIZATION CRITERION        =',1P,E9.2)
 9120 FORMAT(1X,'NEO-CI PROPERTIES WILL BE FOUND FOR ROOT NUMBER',I4)
 9130 FORMAT(1X,'1-PARTICLE DENSITY MATRIX OPTIONS ARE',20I2)
C9140 FORMAT(1X,'PURE SPIN STATE AVERAGED 1- AND 2-PARTICLE DENSITY',
C    *          'MATRIX OPTION=.',1X,'.')
 9150 FORMAT(2(1X,'STATE=',I4,' DM2 WEIGHT=',F10.5,4X,:))
 9155 FORMAT(1X,'STATE-AVERAGED NEO-MCSCF USING IEROOT=',I3,
     *          ' AS ITS STATE-SPECIFIC E')
 9160 FORMAT(/1X,'SYMMETRIES FOR THE',I4,' CORE,',I4,' ACTIVE ARE')
 9200 FORMAT(1X,'THE CURRENT IMPLEMENTAION OF NEO-MCSCF REQUIES THAT'/
     *          'THE ACTIVE NUCLEAR ORBITALS SPAN THE FULL SPACE'/
     *          'I.E., NACT MUST EQUAL THE TOTAL NUMBER OF NUC MOS')
C
      END
C
C*MODULE NEOFCI  *DECK RDCI12_NE
      SUBROUTINE RDCI12_NE(NFT,X2,NCORE,M1,M2,M4,IA,XX,IX,NINTMX,
     *                     M1_NU,MXM2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION X2(M4),IA(MXM2),XX(NINTMX),IX(*)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
C
C     ---- READ TRANSFORMED 1 AND 2 E- INTEGRALS INTO MEMORY ----
C     N.B. THIS ROUTINE IS NOT CORRECT FOR PARALLEL PROCESSING,
C     MESSAGES ARE NEEDED TO MAKE ALL NODES HAVE ALL INTEGRALS
C
      IROW = 0
      DO 110 I=1,MXM2
         IA(I) = IROW
         IROW = IROW+I
  110 CONTINUE
C
C         ONLY THE MASTER HAS THE 1E- INTEGRAL RECORD
C
      CALL SEQREW(NFT)
C
C         READ TRANSFORMED 2E- INTEGRAL FILE IN REVERSE CANONICAL ORDER.
C         TAKE ONLY INTEGRALS WITH ALL 4 INDICES IN THE ACTIVE SPACE.
C
      CALL VCLR(X2,1,M4)
  200 CONTINUE
      CALL PREAD(NFT,XX,IX,NX,NINTMX)
      IF (NX.EQ.0) GO TO 240
      MX = ABS(NX)
      IF (MX.GT.NINTMX) THEN
         IF(MASWRK) WRITE(IW,*) 'CONFUSION WITH INTEGRALS IN -RDCI12-'
         CALL ABRT
         STOP
      END IF
      DO 220 M = 1,MX
         VAL = XX(M)
C
                     NPACK = M
                     IF (LABSIZ .EQ. 2) THEN
*I32                   LABEL1 = IX( 2*NPACK - 1 )
*I32                   LABEL2 = IX( 2*NPACK     )
*I32                   IPACK = ISHFT( LABEL1, -16 )
*I32                   JPACK = IAND( LABEL1, 65535 )
*I32                   KPACK = ISHFT( LABEL2, -16 )
*I32                   LPACK = IAND( LABEL2, 65535 )
*I64                   LABEL = IX(NPACK)
*I64                   IPACK = ISHFT( LABEL, -48 )
*I64                   JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                   KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64                   LPACK = IAND( LABEL, 65535 )
                     ELSE IF (LABSIZ .EQ. 1) THEN
*I32                   LABEL = IX(NPACK)
*I32                   IPACK = ISHFT( LABEL, -24 )
*I32                   JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32                   KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I32                   LPACK = IAND( LABEL, 255 )
*I64                   IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                     LABEL = IX( NPACK/2 )
*I64                     IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64                     JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64                     KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                     LPACK = IAND( LABEL, 255 )
*I64                   ELSE
*I64                     LABEL = IX( (NPACK/2)+1 )
*I64                     IPACK = ISHFT( LABEL, -56 )
*I64                     JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64                     KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64                     LPACK = IAND( ISHFT( LABEL, -32 ), 255 )
*I64                   END IF
                     END IF
                     K = IPACK
                     L = JPACK
                     I = KPACK
                     J = LPACK
C
         I = I-NCORE
         J = J-NCORE
         K = K-NCORE
         L = L-NCORE
         IF(I.LE.0 .OR. I.GT.M1_NU) GO TO 220
         IF(J.LE.0 .OR. J.GT.M1_NU) GO TO 220
         IF(K.LE.0 .OR. K.GT.M1) GO TO 220
         IF(L.LE.0 .OR. L.GT.M1) GO TO 220
C
         IJ = IA(I)+J
         KL = IA(K)+L
CSPW     IJKL = IA(IJ) + KL
         IJKL = (IJ-1)*M2 + KL
         X2(IJKL) = VAL
  220 CONTINUE
      IF (NX .GT. 0) GO TO 200
C
  240 CONTINUE
      CALL SEQREW(NFT)
      RETURN
      END
C
C*MODULE NEOFCI  *DECK CNTDET
      SUBROUTINE CNTDET(NALP,NBLP,NALP_NU,NBLP_NU,IDS,IIS,SOME)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SOME
C
      PARAMETER (MXRT=100)
C
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,S,SZ,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORB,
     *                NA,NB,K,KST,IROOT,IPURES,MAXW1,NITER,MAXP,NCI,
     *                IGPDET,KSTSYM,NFTGCI
      COMMON /DETNEO/ WSTATE_NEO(MXRT),SPINS_NU(MXRT),CRIT_NU,
     *                PRTTOL_NEO,S_NU,SZ_NU,GRPDET_NU,STSYM_NU,
     *                NFLGDM_NEO(MXRT),IWTS_NU(MXRT),NCORSV_NU,
     *                NCOR_NU,NACT_NU,NORB_NU,NA_NU,NB_NU,
     *                K_NEO,KST_NEO,IROOT_NEO,IPURES_NEO,MAXW1_NEO,
     *                NITER_NEO,MAXP_NEO,NCI_NU,IGPDET_NU,KSTSYM_NU
C
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IS,IJKT,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
C
C
      DATA C1/8HC1      /
C
C     --- ELECTRONIC CORE ---
C
      NTOT = NACT + NCORSV
      NTCO = NCORSV
      NCOR = 0
      NSYM = 2**IGPDET
C
C     --- NUCLEAR CORE ---
C
      NTOT_NU = NACT_NU + NCORSV_NU
      NTCO_NU = NCORSV_NU
C
C     --- GET MEMORY NEED TO WORK OUT NUMBER OF DETERMINANTS ---
C         DECIDE NECESSARY DOUBLE/INTEGER WORKING STORAGE -IDS- AND -IIS-
C
      MXNACT = MAX(NACT,NACT_NU)
      MXNTOT = MAX(NTOT,NTOT_NU)
      MXNA   = MAX(NA,NA_NU)
C
      CALL VALFM(LOADFM)
      LIFA  = LOADFM + 1
      LAST  = LIFA   + ((MXNACT+1)*(MXNACT+1))/NWDVAR + 1
      NEED1 = LAST - LOADFM - 1
      CALL GETFM(NEED1)
C
      CALL VALFM(LOADFM)
      IBO = LOADFM + 1
      ICON = IBO + MXNTOT
      ICA  = ICON + MXNA
      ICB  = ICA + NSYM
      KTAB = ICB + NSYM
      IWRK = KTAB + NSYM
      LAST = IWRK + 43
      NEEDT = LAST - LOADFM - 1
      CALL GETFM(NEEDT)
C
C     --- COMPUTE THE TOTAL NUMBER OF ELECTRONIC DETERMINANTS ---
C
      CALL BINOM6(X(LIFA),NACT)
      CALL DAREAD(IDAF,IODA,X(IBO),NTOT,262,1)
      CALL CORTRA(X(IBO),NTOT,NTCO)
      CALL MEMCI(IW,NORB,NCOR,NA,NB,K,MAXP,MAXW1,X(LIFA),
     *           NCI,IDS1,IIS1,NALP,NBLP,IGPDET,KSTSYM,NSYM,
     *           X(IBO),ISST,X(ICON),X(ICA),X(ICB),X(KTAB),X(IWRK))
C
C     --- COMPUTE THE TOTAL NUMBER OF NUCLEAR DETERMINANTS ---
C
      CALL BINOM6(X(LIFA),NACT_NU)
CSPW  NEED NUC MO SYMMETRY LABELS
      CALL DAREAD(IDAF,IODA,X(IBO),NTOT_NU,453,1)
CSPW  CALL PUTINT(X(IBO),NTOT_NU)
      CALL CORTRA(X(IBO),NTOT_NU,NTCO_NU)
      CALL MEMCI(IW,NORB_NU,NCOR_NU,NA_NU,NB_NU,K_NEO,MAXP,MAXW1,
     *           X(LIFA),NCI_NU,IDS2,IIS2,NALP_NU,NBLP_NU,IGPDET,
     *           KSTSYM,NSYM,X(IBO),ISST,X(ICON),X(ICA),X(ICB),
     *           X(KTAB),X(IWRK))
C
      CALL RETFM(NEEDT)
      IDS = MAX(IDS1,IDS2)
      IIS = MAX(IIS1,IIS2)
C
      IF(SOME) THEN
         WRITE(IW,9110) STSYM,GRPDET,SZ,NCI
         IF(GRPDET.NE.C1) GO TO 6
         SVAL = SZ
    5    CONTINUE
            NOCSF = NUMCSF(NA+NB,NACT,SVAL)
            IF(NOCSF.LE.0) GO TO 6
            WRITE(IW,9120) NOCSF,SVAL
            SVAL = SVAL + 1.0D+00
         GO TO 5
      END IF
    6 CONTINUE
C
      IF(SOME) THEN
         WRITE(IW,9115) STSYM,GRPDET,SZ_NU,NCI_NU
         IF(GRPDET.NE.C1) GO TO 10
         SVAL = SZ_NU
    9    CONTINUE
            NOCSF = NUMCSF(NA_NU+NB_NU,NACT_NU,SVAL)
            IF(NOCSF.LE.0) GO TO 10
            WRITE(IW,9120) NOCSF,SVAL
            SVAL = SVAL + 1.0D+00
         GO TO 9
      END IF
   10 CONTINUE
      CALL RETFM(NEED1)
C
      IF(SOME) CALL TIMIT(1)
      RETURN
C
 9110 FORMAT(/1X,'THE NUMBER OF ELECTRONIC DETERMINANTS',
     *           ' HAVING SPACE SYMMETRY ',A2/
     *        1X,'IN POINT GROUP ',A4,' WITH SZ=',F5.1,' IS',I15)
 9115 FORMAT(/1X,'THE NUMBER OF NUCLEAR DETERMINANTS',
     *           ' HAVING SPACE SYMMETRY ',A2/
     *        1X,'IN POINT GROUP ',A4,' WITH SZ=',F5.1,' IS',I15)
 9120 FORMAT(1X,'WHICH INCLUDES',I15,' CSFS WITH S=',F5.1)
      END
C
C*MODULE NEOFCI   *DECK PRIDTN
      SUBROUTINE PRIDTN(ECI,CIVECT,NSTNEO,MIXDET,NACT,NA,NB,
     *                  NACT_NU,NA_NU,NB_NU,PRTTOL,
     *                  IACON,IBCON,IACONN,IBCONN)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION ECI(MIXDET),CIVECT(MIXDET*NSTNEO)
      DIMENSION IFA(0:30,0:30)
      DIMENSION IACON(NA),IBCON(NA),IACONN(NA_NU),IBCONN(NA_NU)
      CHARACTER*30 CONA,CONB,CONNA,CONNB
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
C
      CALL BINOM6(IFA,30)
      NALP = IFA(NACT,NA)
      NBLP = IFA(NACT,NB)
CSPW  NNALP = IFA(NACT_NU,NA_NU)
      NNBLP = IFA(NACT_NU,NB_NU)
C
      IOP=2
      WRITE(IW,9160) PRTTOL
C
      WRITE (IW,*)'FORMAT:FIRST NUCLEAR ALPHA,BETA, THEN ELECTRONIC'
     &            ,' ALPHA,BETA'
      DO I =1,NSTNEO
         WRITE (IW,*)
         WRITE (IW,1100) I,ECI(I)
         WRITE (IW,*)

         CONA(1:30)  = '                              '
         CONB(1:30)  = '                              '
         CONNA(1:30) = '                              '
         CONNB(1:30) = '                              '
C
C        SET UP THE TABLE
C
         IA = (NACT+2)/2 - 2
         IF (IA.LE.0) IA = 1
         CONA(IA:IA+4) = 'ALPHA'
         CONB(IA:IA+4) = 'BETA '
         IAN = (NACT_NU+2)/2 - 2
         IF (IAN.LE.0) IAN = 1
         CONNA(IAN:IAN+4) = 'ALPHA'
         CONNB(IAN:IAN+4) = 'BETA '
         WRITE(IW,'(8A)') CONNA(1:NACT_NU+2),'|',CONNB(1:NACT_NU+2),
     *                    '||',CONA(1:NACT+2),'|',CONB(1:NACT+2),
     *                    '| COEFFICIENT'
         DO 15 II=1,NACT+2
            CONA(II:II) = '-'
15       CONTINUE
         DO 25 II=1,NACT_NU+2
            CONNA(II:II) = '-'
25       CONTINUE
         WRITE(IW,'(8A)') CONNA(1:NACT_NU+2),'|',CONNA(1:NACT_NU+2),
     *                    '||', CONA(1:NACT+2),'|',CONA(1:NACT+2),
     *                    '|------------'
C        NOW PRINT OUT ALL DETERMINANTS
C
         DO 777 K=1,MIXDET
            INN=(I-1)*MIXDET+K
            IF (IOP.EQ.2. AND .ABS(CIVECT(INN)).LT.PRTTOL) GOTO 777
            CALL NDINDX(K,INDA,INDB,INNDA,INNDB,NALP,NBLP,NNBLP)
            DO II=1,NA
               IACON(II) = II
            END DO
            DO II=1,NB
               IBCON(II) = II
            END DO
            DO II=1,NA_NU
               IACONN(II) = II
            END DO
            DO II=1,NB_NU
               IBCONN(II) = II
            END DO
            DO II=1,INDA-1
               CALL ADVANC(IACON,NA,NACT)
            END DO
            DO II=1,INDB-1
               CALL ADVANC(IBCON,NB,NACT)
            END DO
            DO II=1,INNDA-1
               CALL ADVANC(IACONN,NA_NU,NACT_NU)
            END DO
            DO II=1,INNDB-1
               CALL ADVANC(IBCONN,NB_NU,NACT_NU)
            END DO
            CONA(1:30)  = ' 00000000000000000000000000000'
            CONB(1:30)  = ' 00000000000000000000000000000'
            CONNA(1:30) = ' 00000000000000000000000000000'
            CONNB(1:30) = ' 00000000000000000000000000000'
            DO II=1,NA
               CONA(IACON(II)+1:IACON(II)+1) = '1'
            END DO
            DO II=1,NB
               CONB(IBCON(II)+1:IBCON(II)+1) = '1'
            END DO
            DO II=1,NA_NU
               CONNA(IACONN(II)+1:IACONN(II)+1) = '1'
            END DO
            DO II=1,NB_NU
               CONNB(IBCONN(II)+1:IBCONN(II)+1) = '1'
            END DO
C
            CONA(NACT+2:NACT+2) = ' '
            CONB(NACT+2:NACT+2) = ' '
            CONNA(NACT_NU+2:NACT_NU+2) = ' '
            CONNB(NACT_NU+2:NACT_NU+2) = ' '
C
            WRITE(IW,'(8A,F10.7)') CONNA(1:NACT_NU+2),'|',
     *                             CONNB(1:NACT_NU+2),
     *                             '||',CONA(1:NACT+2),'|',
     *                             CONB(1:NACT+2),'|  ',CIVECT(INN)
 777        CONTINUE
      END DO

C
 1100 FORMAT(1X,6HSTATE=,I3,' ENERGY = ',F20.10)
 9160 FORMAT(1X,'PRINTING CI COEFFICIENTS LARGER THAN',F10.6)
      RETURN
      END
C
C*MODULE NEOFCI   *DECK NDINDX
      SUBROUTINE NDINDX(K,INDA,INDB,INNDA,INNDB,NALP,NBLP,NNBLP)
C
      I4=(NNBLP*NALP*NBLP)
      INNDA=(K-1)/I4+1
      IK=MOD(K,I4)
      IF (IK.EQ.0) THEN
         INNDB=NNBLP
         INDA=NALP
         INDB=NBLP
         GOTO 88
      END IF
      I3=(NALP*NBLP)
      INNDB=(IK-1)/I3+1
      IK=MOD(IK,I3)
      IF (IK.EQ.0) THEN
         INDA=NALP
         INDB=NBLP
         GOTO 88
      END IF
      I2=NBLP
      INDA=(IK-1)/I2+1
      INDB=MOD(IK,I2)
      IF (INDB.EQ.0) INDB=NBLP
 88   CONTINUE
C      WRITE (*,*) INNDA,INNDB,INDA,INDB, '  FOR DET # ',K
C
      RETURN
      END
C*MODULE NEOFCI  *DECK RETCIH
      SUBROUTINE RETCIH(CIHAM,SINT1,SINT2,IACON1,IBCON1,IACON2,IBCON2,
     *                 INDEX,NALP,NBLP,NACT,NA,NB,ICOUP,RETCOU,NCOR)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL RETCOU
C
      DIMENSION CIHAM(*),SINT1(*),SINT2(*)
      DIMENSION IACON1(NA),IBCON1(NA),IACON2(NA),IBCON2(NA)
      DIMENSION INDEX(((NACT*(NACT+1))/2),((NACT*(NACT+1))/2))
      DIMENSION ICOUP(*)
C
C     --- SUBROUTINE TO RETURN THE CI HAMILTONIAN MATRIX ---
C
      M1 = NACT
      M2 = (M1*M1+M1)/2
C
C     --- SET UP INDEXING ARRAY ---
C
      DO 5 I=1,M2
         DO 10 J=1,I
            INDEX(I,J) = I*(I-1)/2 + J
            INDEX(J,I) = INDEX(I,J)
   10    CONTINUE
    5 CONTINUE
C
      IHELEM = 0
      DO 20 I=1,NA
         IACON1(I) = I
   20 CONTINUE
C
C     --- OUTER LOOP OVER ALPHA DETERMINANTS ---
C
      ICON1 = 1
      DO 30 II=1,NALP
C
         DO 40 I=1,NB
            IBCON1(I) = I
   40    CONTINUE
C
C        --- LOOP OVER BETA DETERMINANTS ---
C
         DO 50 JJ=1,NBLP
C
            ICON2 = 1
            DO 60 I=1,NA
               IACON2(I) = I
   60       CONTINUE
C
C           --- INNER LOOP OVER ALPHA DETERMINANTS ---
C
            DO 70 KK =1,NALP
C
               DO 80 I=1,NB
                  IBCON2(I) = I
   80          CONTINUE
C
C              --- LOOP OVER BETA DETERMINANTS ---
C
               DO 100 LL=1,NBLP
                  IF(ICON2.GT.ICON1) GOTO 150
C
C                 --- CALCULATE MATRIX ELEMENT ---
C
                  CALL HELEM(SINT1,SINT2,NACT,NA,NB,IACON1,IBCON1,
     *                       IACON2,IBCON2,NCOR,ICON1,ICON2,INDEX,ELEM,
     *                       IONE,IND1,IPERM)
                  IHELEM = IHELEM + 1
                  CIHAM(IHELEM) = ELEM
C
C                 --- STORE COUPLING COEFFICENTS FOR LATER USE ---
C
                  ICP = IND1*IPERM*IONE
                  IF(RETCOU) ICOUP(IHELEM) = ICP
C
                  CALL ADVANC(IBCON2,NB,NACT)
                  ICON2 = ICON2 + 1
  100          CONTINUE
C
               CALL ADVANC(IACON2,NA,NACT)
C
   70       CONTINUE
C
  150       CONTINUE
            CALL ADVANC(IBCON1,NB,NACT)
            ICON1 = ICON1 + 1
C
   50    CONTINUE
C
         CALL ADVANC(IACON1,NA,NACT)
C
  30  CONTINUE
C
      RETURN
      END
C
C*MODULE NEOFCI  *DECK MXHELEM
      SUBROUTINE MXHELEM(SINT2,INDEX,ELEM,
     *                   NA,NB,IA1,IB1,I1,I2,
     *                   NA_NU,NB_NU,IA1_NU,IB1_NU,
     *                   I1_NU,I2_NU,M2,M2_NU,IDXT,IECOUP,
     *                   INCOUP)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION SINT2(*),INDEX(M2,M2_NU),IDXT(*)
      DIMENSION IA1(NA),IB1(NB)
CSPW  DIMENSION IA2(NA),IB2(NB)
      DIMENSION IA1_NU(NA_NU),IB1_NU(NB_NU)
CSPW  DIMENSION IA2_NU(NA_NU),IB2_NU(NB_NU)
      DIMENSION IECOUP(*),INCOUP(*)
C
C     RETURNS THE MATRIX ELEMENT < IK | H(MIX) |J L > WHERE K, L ARE
C     DETERMINANTS.  ALPHA AND BETA OCCUPATIONS ARE STORED IN
C     IA1, IB1 FOR K AND IA2, IB2 FOR L.
C
C     THIS CODE PRESENTLY USES SAVED COUPLING COEFFICENTS. THE NEXT
C     IMPEMENTATION WILL HAVE THE ABILIY TO FORM THEM ON THE FLY. SPW.
C
      ELEM = 0.0D+00
C
      IF(I1.GT.I2) THEN
         IE  = I1*(I1-1)/2 + I2
      ELSE
         IE  = I2*(I2-1)/2 + I1
      END IF
      IF(I1_NU.GT.I2_NU) THEN
         IN = I1_NU*(I1_NU-1)/2 + I2_NU
      ELSE
         IN = I2_NU*(I2_NU-1)/2 + I1_NU
      END IF
C
C     --- IF ELEC DETS THE SAME AND NUC DETS THE SAME ---
C
      IF (I1.EQ.I2 .AND. I1_NU.EQ.I2_NU) THEN
C
C        --- ALPHA ELECS - ALPHA NUC ---
         DO 10 I=1,NA
            IA = IDXT(IA1(I))
            DO 20 J=1,NA_NU
               JA = IDXT(IA1_NU(J))
               IX = INDEX(IA,JA)
               ELEM = ELEM - SINT2(IX)
   20       CONTINUE
   10    CONTINUE
C
C        --- BETA ELECS - ALPHA NUC ---
         DO 30 I=1,NB
            IB = IDXT(IB1(I))
            DO 40 J=1,NA_NU
               JA = IDXT(IA1_NU(J))
               IX = INDEX(IB,JA)
               ELEM = ELEM - SINT2(IX)
   40       CONTINUE
   30    CONTINUE
C
C        --- ALPHA ELECS - BETA NUC ---
         DO 50 I=1,NA
            IA = IDXT(IA1(I))
            DO 60 J=1,NB_NU
               JB = IDXT(IB1_NU(J))
               IX = INDEX(IA,JB)
               ELEM = ELEM - SINT2(IX)
   60       CONTINUE
   50    CONTINUE
C
C        --- BETA ELECS - BETA NUC ---
         DO 70 I=1,NB
            IB = IDXT(IB1(I))
            DO 80 J=1,NB_NU
               JB = IDXT(IB1_NU(J))
               IX = INDEX(IB,JB)
               ELEM = ELEM - SINT2(IX)
   80       CONTINUE
   70    CONTINUE
C
C
C     --- IF ELECTRONIC DETERMINANTS ARE SAME ---
C
      ELSE IF (I1.EQ.I2) THEN
C
         NUIDX = INCOUP(IN)
         IF(NUIDX.EQ.0) GOTO 120
         JA = ABS(NUIDX)
         SIGN = DBLE(NUIDX/JA)
C
C        --- ALPHA ELECS ---
         DO 100 I=1,NA
            IA = IDXT(IA1(I))
            IX = INDEX(IA,JA)
            ELEM = ELEM - SIGN*SINT2(IX)
  100    CONTINUE
C
C        --- BETA ELECS ---
         DO 110 I=1,NB
            IB = IDXT(IB1(I))
            IX = INDEX(IB,JA)
            ELEM = ELEM - SIGN*SINT2(IX)
  110    CONTINUE
C
  120    CONTINUE
C
C     --- IF NUCLEAR DETERMINANTS ARE SAME ---
C
      ELSE IF (I1_NU.EQ.I2_NU) THEN
C
         IDXEL = IECOUP(IE)
         IF(IDXEL.EQ.0) GOTO 160
         IA = ABS(IDXEL)
         SIGN = DBLE(IDXEL/IA)
C
C        --- ALPHA NUC ---
         DO 140 J=1,NA_NU
            JA = IDXT(IA1_NU(J))
            IX = INDEX(IA,JA)
            ELEM = ELEM - SIGN*SINT2(IX)
  140    CONTINUE
C
C        ---  BETA NUC ---
         DO 150 J=1,NB_NU
            JB = IDXT(IB1_NU(J))
            IX = INDEX(IA,JB)
            ELEM = ELEM - SIGN*SINT2(IX)
  150    CONTINUE
C
  160    CONTINUE
C
C     --- BOTH ELEC AND NUCLEAR DETERMINANTS DIFFERENT ---
C
      ELSE
C
         NUIDX = INCOUP(IN)
         IF(NUIDX.EQ.0) GOTO 200
         JA = ABS(NUIDX)
         SIGN1 = DBLE(NUIDX/JA)
C
         IDXEL = IECOUP(IE)
         IF(IDXEL.EQ.0) GOTO 200
         IA = ABS(IDXEL)
         SIGN2 = DBLE(IDXEL/IA)
C
         IX = INDEX(IA,JA)
         ELEM = ELEM - SIGN1*SIGN2*SINT2(IX)
C
  200    CONTINUE
C
      ENDIF
C
      RETURN
      END
C
C*MODULE NEOFCI  *DECK DETFCI_RETMAT
      SUBROUTINE DETFCI_RETMAT(CLABEL,SINT1,SINT2,IA,
     *           XX,IXX,IACON1,IBCON1,IACON2,
     *           IBCON2,IACON1_NU,IBCON1_NU,IACON2_NU,
     *           IBCON2_NU,INDEX,HNMO,CIEHAM,IECOUP,CINHAM,
     *           INCOUP,CINEH,NALP,NBLP,NALP_NU,NBLP_NU,
     *           M1,M2,M4,M1_NU,M2_NU,M4_NU,M5,MXM1,MXM2,
     *           MAX2IN,NINTMX,MXNA,NCORSV,NACT,
     *           NA,NB,NCI,NACT_NU,NA_NU,NB_NU,NCI_NU)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL CLABEL
C
      COMMON /IOFILE/ IR,IW,IP,IS,IJKT,IDAF,NAV,IODA(950)
C
      DIMENSION SINT1(MXM2),SINT2(MAX2IN)
      DIMENSION IA(*),XX(NINTMX),IXX(NINTMX)
      DIMENSION IACON1(MXNA),IBCON1(MXNA),IACON2(MXNA),IBCON2(MXNA)
      DIMENSION IACON1_NU(MXNA),IBCON1_NU(MXNA)
      DIMENSION IACON2_NU(MXNA),IBCON2_NU(MXNA)
      DIMENSION INDEX(*),HNMO(*)
      DIMENSION CIEHAM((NCI*NCI+NCI)/2)
      DIMENSION IECOUP((NCI*NCI+NCI)/2)
      DIMENSION CINHAM((NCI_NU*NCI_NU+NCI_NU)/2)
      DIMENSION INCOUP((NCI_NU*NCI_NU+NCI_NU)/2)
      DIMENSION CINEH(*)
C
C     --- FORM MIXED NUCLEAR-ELECTRONIC CI (NEO-CI) MATRIX ELEMENTS ---
C
C     --- CORE CONTRIBUTION TO THE ENERGY IS OBTAINED FROM -ECORE-, ---
C         AND FROM MODIFICATIONS TO THE TRANSFORMED 1E- INTEGRALS.
C         THIS EFFECTIVELY REMOVES CORE ORBITALS FROM THE COMPUTATION.
C
CSPW  USE OF SYMMETRY NOT YET CODED
C
C     -- OBTAIN 1 AND 2 E- TRANSFORMED INTEGRALS OVER ACTIVE ORBITALS --
C     CALLING ARGUMENT -CLABEL- GOVERNS WHETHER TRANSFORMED INTEGRALS
C     ON FILE -IJKT- INCLUDE THE CORE ORBITALS OR NOT.  IT IS ASSUMED
C     THAT NO CORE INTEGRALS ARE IN -IJKT-, SO THIS VARIABLE TELLS IF
C     THE ACTIVE ORBITALS START FROM 1,2,3... OR NCORSV+1,NCORSV+2,...
C
      NCORE = 0
      IF(CLABEL) NCORE = NCORSV
      CALL RDCI12(.FALSE.,IJKT,SINT1,SINT2,NCORE,M1,M2,M4,
     *            IA,XX,IXX,NINTMX)
C
C     --- COMPUTE ELECTRONIC CI HAMILTONIAN MATRIX ELEMENTS HE(IJ) ---
C
      CALL RETCIH(CIEHAM,SINT1,SINT2,IACON1,IBCON1,IACON2,IBCON2,
     *            INDEX,NALP,NBLP,NACT,NA,NB,IECOUP,.TRUE.,0)
C
C
      IF(NA_NU.EQ.1 .AND. NB_NU.EQ.0) THEN
C
C        --- IF ONLY ONE QUAUNTUM NUCLEUS NUCLEAR CI HAMILTONIAN IS
C            SIMPLY EQUAL TO THE 1-NUCLEAR INTEGRALS IN THE MO ---
C
         CALL DAREAD(IDAF,IODA,SINT1,M2_NU,459,0)
         CALL DCOPY(M2_NU,SINT1,1,CINHAM,1)
         IF(NCORSV.GT.0) THEN
            CALL DAREAD(IDAF,IODA,HNMO,M2_NU,460,0)
            CALL VADD(CINHAM,1,HNMO,1,CINHAM,1,M2_NU)
         END IF
C
C        --- SET UP 'COUPLING COEFFICIENTS' FOR USE IN MIXED HAM ---
C        (IN THIS SPECIAL CASE JUST INDEXES OF THE OFF DIAGONALS)
C
         IJ = 0
         DO I=1,NCI_NU
            DO J=1,I
                IJ = IJ + 1
                INCOUP(IJ) = IJ
                IF(I.EQ.J) INCOUP(IJ) = 0
            END DO
         END DO
C
      ELSE
C
C     --- OBTAIN 1 AND 2 NUC TRANSFORMED INTEGRALS OVER ACTIVE
C         ORBITALS --
C
         NCORN = 0
CSPW     IF(CLABEL) NCORN = NCORSV_NU
         NNIJKT = 69
         CALL RDCI12(.FALSE.,NNIJKT,SINT1,SINT2,NCORN,M1_NU,M2_NU,
     *               M4_NU,IA,XX,IXX,NINTMX)
C
C        --- 1-NUC INTEGRALS MODIFIED DUE TO ELECTRONIC CORE ---
C
         IF(NCORSV.GT.0) THEN
            CALL DAREAD(IDAF,IODA,HNMO,M2_NU,460,0)
            CALL VADD(SINT1,1,HNMO,1,SINT1,1,M2_NU)
         END IF
C
C        --- COMPUTE NUCLEAR CI HAMILTONIAN MATRIX ELEMENTS HN(KL) ---
C
         CALL RETCIH(CINHAM,SINT1,SINT2,IACON1,IBCON1,
     *               IACON2,IBCON2,INDEX,NALP_NU,NBLP_NU,NACT_NU,
     *               NA_NU,NB_NU,INCOUP,.TRUE.,0)
C
      END IF
C
C     --- GET ELECTRON-NUCLEAR INTEGRALS OVER ACTIVE ORBITALS ---
C
CSPW  NCORE = 0
CSPW  IF(CLABEL) NCORE = NCORSV
      NEIJKT = 71
      NCORE = 0
      CALL RDCI12_NE(NEIJKT,SINT2,NCORE,M1,M2,M5,IA,
     *               XX,IXX,NINTMX,M1_NU,MXM2)
C
C     --- NOW FORM NUCLEAR-ELECTRONIC CI (NEO-CI) HAMILTONIAN ELEMENTS ---
C         ACCORDING TO:
C         H(IK),(JL) = HE(IJ)<K|L> + HN(KL)<I|J> - (CROSS TERM)
C
      CALL NEOCIH(CIEHAM,IECOUP,CINHAM,INCOUP,CINEH,SINT2,INDEX,
     *            IACON1,IBCON1,IACON2,IBCON2,
     *            NALP,NBLP,NACT,NA,NB,
     *            IACON1_NU,IBCON1_NU,IACON2_NU,IBCON2_NU,
     *            NALP_NU,NBLP_NU,NACT_NU,NA_NU,NB_NU,M2,
     *            M2_NU,IA,MXM1)
C
      RETURN
C
      END
C*MODULE NEOFCI  *DECK NEOCIH
      SUBROUTINE NEOCIH(CIEHAM,IECOUP,CINHAM,INCOUP,CINEH,SINT2,INDEX,
     *                  IA1,IB1,IA2,IB2,NALP,NBLP,NACT,NA,NB,
     *                  IA1_NU,IB1_NU,IA2_NU,IB2_NU,NALP_NU,NBLP_NU,
     *                  NACT_NU,NA_NU,NB_NU,
     *                  M2,M2_NU,IA,IALIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION CIEHAM(*),IECOUP(*),CINHAM(*),INCOUP(*)
      DIMENSION CINEH(*),SINT2(*),INDEX(M2,M2_NU),IA(*)
      DIMENSION IA1(NA),IB1(NA),IA2(NA),IB2(NA)
      DIMENSION IA1_NU(NA_NU),IB1_NU(NA_NU),IA2_NU(NA_NU),IB2_NU(NA_NU)
C
C     --- ROUTINE TO RETURN THE NEO-CI HAMILTONIAN MATRIX ---
C
C     ON ENTRY: CIEHAM  - ELECTRONIC CI HAMILTONIAN MATRIX
C               IECOUP  - 1-E COUPLING COEFFICENTS
C               CINHAM  - NUCLEAR ELECTRONIC CI HAMILTONIAN
C               INCOUP  - 1-N COUPLING COEFFICENTS
C               SINT2   - ELECTRONIC-NUCLEAR MO INTEGRALS
C               INDEX   - ARRAY TO BE FILLED WITH INTEGRAL INDEXING
C
C     ON EXIT:  CINEH  - NUCLEAR-ELECTRONIC CI HAMILTONIAN MATRIX
C
C
C     --- SET UP INDEXING ARRAY (IJ=ACTIVE NUC; KL=ACTIVE ELECS) ---
C
      DO 5 KL=1,M2
         DO 10 IJ=1,M2_NU
            INDEX(KL,IJ) = (IJ-1)*M2 + KL
   10    CONTINUE
    5 CONTINUE
C
      DO II=1,IALIM
         IA(II) = (II*II+II)/2
      END DO
C
      IHELEM = 0
C
      DO 20 I=1,NA_NU
         IA1_NU(I) = I
   20 CONTINUE
C
C     --- OUTER LOOP OVER ALPHA NUCLEAR DETERMINANTS ---
C
      IMIX1 = 1
      I1_NU = 0
      DO 30 IIA=1,NALP_NU
C
         DO 40 I=1,NB_NU
            IB1_NU(I) = I
   40    CONTINUE
C
C        --- OUTER LOOP OVER BETA NUCLEAR DETERMINANTS ---
C
         DO 50 IIB=1,NBLP_NU
            I1_NU = I1_NU + 1
C
            I1 = 0
            DO 60 I=1,NA
               IA1(I) = I
   60       CONTINUE
C
C           --- OUTER LOOP OVER ALPHA ELECTRONIC DETERMINANTS ---
C
            DO 70 KKA=1,NALP
C
               DO 80 I=1,NB
                  IB1(I) = I
   80          CONTINUE
C
C              --- OUTER LOOP OVER BETA ELECTRONIC DETERMINANTS ---
C
               DO 100 KKB=1,NBLP
                  I1 = I1 + 1
C
                  I2_NU = 0
                  DO 110 I=1,NA_NU
                     IA2_NU(I) = I
  110             CONTINUE
C
C                 --- INNER LOOP OVER ALPHA NUCLEAR DETERMINANTS ---
C
                  IMIX2 = 1
                  DO 120 JJA=1,NALP_NU
C
                     DO 130 I=1,NB_NU
                        IB2_NU(I) = I
  130                CONTINUE
C
C                    --- LOOP OVER BETA NUCLEAR DETERMINANTS ---
C
                     DO 140 JJB=1,NBLP_NU
                        I2_NU = I2_NU + 1
C
                        I2 = 0
                        DO 150 I=1,NA
                           IA2(I) = I
  150                   CONTINUE
C
C                       --- INNER LOOP OVER ALPHA ELECTRONIC DETERMINANTS ---
C
                        DO 160 LLA=1,NALP
C
                           DO 170 I=1,NB
                              IB2(I) = I
  170                      CONTINUE
C
C                          -- INNER LOOP OVER BETA ELECTRONIC DETERMINANTS --
C
                           DO 180 LLB=1,NBLP
                              IF(IMIX2.GT.IMIX1) GOTO 350
                              I2 = I2 + 1
C
C                             --- GET NUC-ELEC CROSS TERM ELEMENTS ---
C
                              IHELEM = IHELEM + 1
                              CALL MXHELEM(SINT2,INDEX,ELEM,
     *                             NA,NB,IA1,IB1,I1,I2,
     *                             NA_NU,NB_NU,IA1_NU,IB1_NU,
     *                             I1_NU,I2_NU,
     *                             M2,M2_NU,IA,IECOUP,INCOUP)
C
C
C                             --- ADD NUC CI MATRIX ELEMENTS ---
C
                              IF(I1.EQ.I2) THEN
                                 IF(I1_NU.GT.I2_NU) THEN
                                    IN = I1_NU*(I1_NU-1)/2 + I2_NU
                                 ELSE
                                    IN = I2_NU*(I2_NU-1)/2 + I1_NU
                                 END IF
                                 ELEM = ELEM + CINHAM(IN)
                              END IF
C
C                             --- ADD ELEC CI MATRIX ELEMENTS ---
C
                              IF(I1_NU.EQ.I2_NU) THEN
                                 IF(I1.GT.I2) THEN
                                    IE  = I1*(I1-1)/2 + I2
                                 ELSE
                                    IE  = I2*(I2-1)/2 + I1
                                 END IF
                                 ELEM = ELEM + CIEHAM(IE)
                              END IF
C
C                             --- ELEMENT COMPLETE PUT IN MATRIX ---
C
                              CINEH(IHELEM) = ELEM
C
                              CALL ADVANC(IB2,NB,NACT)
                              IMIX2 = IMIX2 + 1
  180                      CONTINUE
C
                           CALL ADVANC(IA2,NA,NACT)
C
  160                   CONTINUE
C
                        CALL ADVANC(IB2_NU,NB_NU,NACT_NU)
C
  140                CONTINUE
C
                     CALL ADVANC(IA2_NU,NA_NU,NACT_NU)
C
  120             CONTINUE
C
  350             CONTINUE
                  CALL ADVANC(IB1,NB,NACT)
                  IMIX1 = IMIX1 + 1
C
  100          CONTINUE
C
               CALL ADVANC(IA1,NA,NACT)
C
   70       CONTINUE
C
            CALL ADVANC(IB1_NU,NB_NU,NACT_NU)
C
   50    CONTINUE
         CALL ADVANC(IA1_NU,NA_NU,NACT_NU)
C
   30 CONTINUE
C
      RETURN
      END
C
C*MODULE NEOFCI  *DECK EDETIN
      SUBROUTINE EDETIN(NPRINT,GPNAME)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SOME,PURES,GOPARR,DSKWRK,MASWRK,ABEL,WTSOK,ANALYS,CLOBBR
C
      PARAMETER (MXATM=2000, MXRT=100, MXAO=8192, MXSH=5000)
C
      COMMON /DETPAR/ ICLBBR,ANALYS
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,S,SZ,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORB,
     *                NA,NB,K,KST,IROOT,IPURES,MAXW1,NITER,MAXP,NCI,
     *                IGPDET,KSTSYM,NFTGCI
      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 /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,MA,MB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PSILVL/ IPSI,ISKPRP
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMMOL/ GROUP,COMPLEX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
C
      DIMENSION FANT(8),LFANT(8),GANT(27),LGANT(8)
C
      PARAMETER (NNAM=20)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA QNAM/8HGROUP   ,8HISTSYM  ,
     *          8HNCORE   ,8HNACT    ,8HNELS    ,8HSZ      ,
     *          8HNSTATE  ,8HNSTGSS  ,8HNHGSS   ,8HMXXPAN  ,
     *          8HITERMX  ,8HCVGTOL  ,8HPRTTOL  ,
     *          8HIROOT   ,8HNFLGDM  ,8HPURES   ,8HWSTATE  ,
     *          8HWTSOK   ,8HANALYS  ,8HCLOBBR  /
      DATA KQNAM/5,1,   1,1,1,3,  1,1,1,1,   1,3,3,  1,-1,0,-3,   0,0,0/
C
      DATA RNONE/8HNONE    /
      DATA HESS/8HHESSIAN /
      DATA FANT/8HC1      ,8HCI      ,8HCS      ,8HC2      ,
     *          8HD2      ,8HC2V     ,8HC2H     ,8HD2H     /
      DATA LFANT/1,1,1,1,2,2,2,3/
      DATA GANT/8HA       ,8HAG      ,8HAU      ,8HA'      ,
     *          8HA"      ,8HA       ,8HB       ,8HA       ,
     *          8HB1      ,8HB2      ,8HB3      ,8HA1      ,
     *          8HA2      ,8HB1      ,8HB2      ,8HAG      ,
     *          8HBG      ,8HBU      ,8HAU      ,8HAG      ,
     *          8HB1G     ,8HB2G     ,8HB3G     ,8HAU      ,
     *          8HB1U     ,8HB2U     ,8HB3U     /
      DATA LGANT/0,1,3,5,7,11,15,19/
C
      SOME = MASWRK  .AND.  NPRINT.NE.-5  .AND.  NPRINT.NE.-23
C
      IF(SOME) WRITE(IW,9000)
C
C          SET UP INPUT TO SPECIFY THE CI SPACE
C          SELECT POINT GROUP AND ELECTRON/ORBITAL COUNTS
C
C          THIS ROUTINE IS DETINP WITH VERY MINOR MODIFICATIONS - SPW
C
      GRPDET = FANT(1)
      IF(IGROUP.EQ.1)                GRPDET = FANT(1)
      IF(IGROUP.EQ.3)                GRPDET = FANT(2)
      IF(IGROUP.EQ.2)                GRPDET = FANT(3)
      IF(IGROUP.EQ.4.AND.NAXIS.EQ.2) GRPDET = FANT(4)
      IF(IGROUP.EQ.8.AND.NAXIS.EQ.2) GRPDET = FANT(5)
      IF(IGROUP.EQ.7.AND.NAXIS.EQ.2) GRPDET = FANT(6)
      IF(IGROUP.EQ.6.AND.NAXIS.EQ.2) GRPDET = FANT(7)
      IF(IGROUP.EQ.9.AND.NAXIS.EQ.2) GRPDET = FANT(8)
      IF(NT.EQ.1) GRPDET=FANT(1)
      IF(RUNTYP.EQ.HESS  .AND.  NHLEVL.GT.0) GRPDET=FANT(1)
      KSTSYM = 1
C
      NCORE  = 0
      NACT   = 0
      NELS   = 0
      SZ     = (MUL-1)/TWO
C
C          SET UP INPUT TO CONTROL THE DIAGONALIZATION
C
      NSTATE = 1
      NSTGSS = 1
      NHGSS  = 300
      MXXPAN = 10
      ITERMX = 100
      CVGTOL = 1.0D-05
      CALL DERCHK(NDER)
      IF(NDER.GT.0)      CVGTOL=1.0D-06
      IF(MPLEVL.GT.0)    CVGTOL=1.0D-06
      IF(CITYP.NE.RNONE) CVGTOL=1.0D-06
      IF(NFG.NE.0)       CVGTOL=1.0D-06
      PRTTOL = 0.05D+00
C
C          SET UP INPUT TO CONTROL THE FIRST ORDER DENSITY COMPUTATION
C          MCSCF USES IROOT TO SELECT A SPECIFIC STATE'S E, IF AVERAGING
C
      IROOT=0
      IF(IPSI.EQ.1) IROOT=1
      KQNAM(15)=MXRT*10 + 1
      DO 5 I=1,MXRT
         NFLGDM(I) = 0
    5 CONTINUE
      NFLGDM(1)=1
C
C          SET UP INPUT TO CONTROL THE SECOND ORDER DENSITY COMPUTATION
C
      PURES = .TRUE.
      KQNAM(17)=MXRT*10 + 3
      CALL VCLR(WSTATE,1,MXRT)
      WSTATE(1) = ONE
      WTSOK  = .FALSE.
      ANALYS = .FALSE.
      CLOBBR = .FALSE.
C
      CALL NAMEIO(IR,JRET,GPNAME,NNAM,QNAM,KQNAM,
     *            GRPDET,KSTSYM,NCORE,NACT,NELS,SZ,NSTATE,NSTGSS,
     *            NHGSS,MXXPAN,ITERMX,CVGTOL,PRTTOL,IROOT,NFLGDM,
     *            PURES,WSTATE,WTSOK,ANALYS,CLOBBR,
     *            0,0,0,0,
     *    0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,
     *    0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0)
C
C        FOR FMO-ROHF RUNS USING MCSCF THE ACTIVE SPACE IS
C        AUTOMATICALLY SET, NO INPUT GROUP NEEDED HERE.
C
      IF(JRET.EQ.2.AND.(NFG.EQ.0.OR.MONCOR.GE.0)) THEN
         IF(MASWRK) WRITE(IW,9010) GPNAME
         CALL ABRT
      END IF
C
C     --- NEO CANNOT USE SYMMETRY YET - OVERRIDE ANY THE USER HAS INPUT ---
      GRPDET = FANT(1)
      KSTSYM = 1
C
      IGPDET = -1
C
C        NUMERICAL HESSIANS FOR MCSCF MAKE ATOMIC DISPLACEMENTS INTO C1
C
      IF ((RUNTYP.EQ.HESS.AND.NHLEVL.GT.0)  .OR.  NT.EQ.1) THEN
         GRPDET=FANT(1)
         KSTSYM=1
      END IF
C
C        THE INPUT FOR C2H IS SUPPOSED TO BE IDENTICAL TO THE GUGA
C        ORDER, NAMELY 1,2,3,4=AG,BU,BG,AU, BUT THE CI CODE WANTS
C        THE ORDER OF  1,2,3,4=AG,BG,BU,AU.  SEE ALSO GAJASW ROUTINE.
C
      IF (GRPDET.EQ.FANT(7)) THEN
         MODI = KSTSYM
         IF(KSTSYM.EQ.2) MODI=3
         IF(KSTSYM.EQ.3) MODI=2
         KSTSYM=MODI
      END IF
      DO I=1,8
         IF (GRPDET.EQ.FANT(I)) THEN
            IGPDET=LFANT(I)
            STSYM = GANT(LGANT(I)+KSTSYM)
         ENDIF
      END DO
      IF (IGPDET.EQ.-1) THEN
         IF(MASWRK) WRITE(IW,*) '$DET POINT GROUP IS UNRECOGNIZED!'
         CALL ABRT
      ENDIF
      IF (GRPDET.EQ.FANT(1).AND.KSTSYM.GT.1) THEN
         IF(MASWRK) WRITE(IW,*) '$DET STATE SYMM IS NOT CORRECT IRREP'
         CALL ABRT
      ENDIF
      IF (KSTSYM.GT.(2**IGPDET)) THEN
         IF(MASWRK) WRITE(IW,*)
     *       '$DET STATE SYMMETRY IS TOO LARGE FOR THIS GROUP'
         CALL ABRT
      ENDIF
C
C READ MO SYMMETRIES AND WRITE TO DIRECT ACCESS FILE,
C BUT DON'T DO IF THIS IF WE ARE PEEKING AT THE $DET INPUT
C AT THE MOST EARLY STAGES OF AN INITIAL MCSCF RUN.
C
      IF (NPRINT.EQ.-23) GOTO 1314
C
      L0 = NQMT
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 =  L1*L1
C
      CALL VALFM(LOADFM)
      LMOLAB = LOADFM + 1
      LMOIRP = LMOLAB + L1
      LVEC   = LMOIRP + L1
      LS     = LVEC   + L3
      LQ     = LS     + L2
      LWRK   = LQ     + L3
      LMODEG = LWRK   + L1
      LAST   = LMODEG + L1
      NEEDD   = LAST - LOADFM - 1
      CALL GETFM(NEEDD)
C
      IF(GRPDET.EQ.FANT(1)) THEN
         CALL C1DET(X(LMOIRP),X(LMOLAB),L0)
      ELSE
         CALL DAREAD(IDAF,IODA,X(LVEC),L3,15,0)
         CALL DAREAD(IDAF,IODA,X(LS),L2,12,0)
         CALL DAREAD(IDAF,IODA,X(LQ),L3,45,0)
         CALL TRFSYM(X(LMOLAB),X(LMOIRP),X(LMODEG),X(LQ),X(LS),X(LVEC),
     *               X(LWRK),IA,L0,L1,L0,L1)
      END IF
C
C  CHANGE ORBITAL SYMMETRY LABELS FROM GAMESS TO JAKAL VALUES
C
      CALL GAJASW(X(LMOIRP),NUM,GRPDET)
C
      CALL DAWRIT(IDAF,IODA,X(LMOIRP),L1,262,1)
C
C     1.  SET NCORSV,NCOR,NACT,NORB,NA,NB FOR DETERMINANT SPECIFICATION
C         CHECK INPUT, AND COPY INTO INTERNALLY USED VARIABLE NAMES.
C         NCOR WILL BE SET TO ZERO TO DROP CORES, SO NCORSV SAVES THIS.
C
 1314 CONTINUE
C
C     ADD SECOND RHF MONOMER AS CORE FOR FMO RUNS OF MCSCF DIMERS
C     EMULATE ROHF WITH AN APPROPRIATE MCSCF IF MONCOR IS NEGATIVE.
C
      IF(NFG.NE.0) THEN
         IF(MONCOR.GE.0) THEN
            IF(IFMOSTP.GE.3) NCORE=NCORE+MONCOR
         ELSE
            IF(MUL.EQ.1) MUL=3
            NACT=MUL-1
            NELS=NACT
            NCORE=(NE-NELS)/2
            SZ     = (MUL-1)/TWO
C           WRITE(6,*) 'FMO-ROHF RESET CAS SPACE TO',NACT,NCORE,NELS,SZ
C           CALL ABRT
         ENDIF
      ENDIF
C
      IF(NPRINT.NE.-23  .AND.
     *       (NCORE.LT.0  .OR.  NACT.LE.0  .OR.  NELS.LE.0)) THEN
         IF(MASWRK) WRITE(IW,9020) GPNAME,NCORE,NACT,NELS
         CALL ABRT
      END IF
      NCORSV = NCORE
      NCOR   = NCORE
      NORB   = NCORE + NACT
      NHIGH = INT(SZ+SZ+0.0001D+00)
      NB = (NELS-NHIGH)/2
      NA = NB+NHIGH
      MA = NA+NCORSV
      MB = NB+NCORSV
      NELTOT = 2*NCOR+NA+NB
      NERR=0
      IF(NELTOT.NE.NE)  NERR=1
      IF(NELS.NE.NA+NB) NERR=1
      IF(NA.LT.NB)      NERR=1
      IF(NA.LE.0)       NERR=1
      IF(NB.LT.0)       NERR=1
      IF(NPRINT.NE.-23  .AND.  NERR.GT.0) THEN
         IF(MASWRK) WRITE(IW,9030) NCORE,NELS,SZ,ICH,MUL
         CALL ABRT
      END IF
      S = (MUL-1)/TWO
C
C        2. SET K,KST,MAXW1,NITER,MAXP,CRIT,PRTTOL FOR DIAGONALIZATION
C
      K     = NSTATE
      KST   = MAX(NSTGSS,K)
      MAXP  = MAX(MXXPAN,2*KST)
      MAXW1 = NHGSS
      NITER = ITERMX
      CRIT  = CVGTOL
      IF(NPRINT.NE.-23  .AND.  K.GT.MXRT) THEN
         IF(MASWRK) WRITE(IW,9035) K,MXRT
         CALL ABRT
      END IF
                 ICLBBR=0
      IF(CLOBBR) ICLBBR=1
C
C        3. SETUP FOR 1ST ORDER DENSITY COMPUTATION
C
      IF(IROOT.GT.NSTATE) THEN
         IF(MASWRK) WRITE(IW,9036) IROOT,NSTATE
         CALL ABRT
      END IF
      IF(IROOT.GT.MXRT) THEN
         IF(MASWRK) WRITE(IW,9037) IROOT,MXRT
         CALL ABRT
      END IF
      IF(NFLGDM(IROOT).EQ.0) NFLGDM(IROOT)=1
C
C        4. SETUP FOR STATE-AVERAGING 1ST AND 2ND ORDER DENSITIES.
C        -IWTS- INDEXES THE NON-ZERO ELEMENTS OF -WSTATE-
C
      IPURES=0
      IF(PURES) IPURES=1
      MXSTAT=0
      WSUM = ZERO
      DO 15 I=1,MXRT
         IF(WSTATE(I).GT.ZERO) THEN
            IF(I.LE.NSTATE) THEN
               WSUM = WSUM + WSTATE(I)
               MXSTAT = MXSTAT+1
               IWTS(MXSTAT) = I
            ELSE
               IF(MASWRK) WRITE(IW,9040) NSTATE
               CALL ABRT
            END IF
         END IF
         IF(WSTATE(I).LT.ZERO) THEN
            IF(MASWRK) WRITE(IW,9050)
            CALL ABRT
         END IF
   15 CONTINUE
      SCALE = ONE/WSUM
      CALL DSCAL(MXRT,SCALE,WSTATE,1)
C
      IF(IROOT.GT.MXSTAT) THEN
         IF(MASWRK) WRITE(IW,9055) IROOT,MXSTAT
         CALL ABRT
      END IF
C
C        IF RUNNING SILENTLY, RETURN WITHOUT PRINTING ANYTHING
C
      IF(NPRINT.EQ.-23) RETURN
C
C     THE COMPUTATION OF THE ENERGY GRADIENT REQUIRES A SYMMETRIC
C     LAGRANGIAN (SINCE GAMESS DOES NOT AT PRESENT DO THE CPHF
C     PROBLEM TO OBTAIN ORBITAL PERTURBATIONS).  ALTHOUGH A STATE
C     AVERAGED MCSCF RUN RESULTS IN A SYMMETRIC "AVERAGE LAGRANGIAN",
C     THE LAGRANGIAN FOR ANY SINGLE STATE USING THE AVERAGED ORBITALS
C     WILL *NOT* BE SYMMETRIC.  HENCE ANY JOB THAT ATTEMPTS TO STATE
C     AVERAGE AND ALSO DO A NUCLEAR GRADIENT SHOULD BE FLUSHED.
C     THERE IS A SECRET INPUT KEYWORD TO AVOID BEING FLUSHED, WHICH
C     MAY BE REASONABLE IF THE STATES BEING AVERAGED ARE DEGENERATE.
C
      CALL DERCHK(MAXDER)
      IF(MXSTAT.GT.1  .AND.  MAXDER.GT.0) THEN
         IF(MAXDER.EQ.1  .AND.  NGLEVL.NE.1) THEN
            IF(WTSOK) THEN
               IF (SOME) WRITE(IW,9061) RUNTYP,MXSTAT
            ELSE
               IF (SOME) WRITE(IW,9062) RUNTYP,MXSTAT
               CALL ABRT
            END IF
         END IF
         IF(MAXDER.EQ.2  .AND.  NHLEVL.NE.2) THEN
            IF (SOME) WRITE(IW,9063) RUNTYP,MXSTAT
            CALL ABRT
         END IF
         IF(IROOT.EQ.0   .AND.  NGLEVL+NHLEVL.GT.0) THEN
            IF(SOME) WRITE(IW,9064)
            CALL ABRT
         END IF
      END IF
C
      IF(SOME) THEN
         WRITE(IW,9100) GRPDET,STSYM,NCOR,NACT,
     *                  NA+NCOR,NA,NB+NCOR,NB,NORB
      END IF
C
      IF(SOME) THEN
         WRITE(IW,9160) NCOR,NACT
         CALL MOSYPR(X(LMOLAB),NCOR,NACT)
      ENDIF
C
      CALL RETFM(NEEDD)
      RETURN
C
C
 9000 FORMAT(/5X,23(1H-)/
     *       5X,'ELECTRONIC ACTIVE SPACE'/
     *       5X,23(1H-))
 9010 FORMAT(/1X,'**** ERROR, THIS RUN REQUIRES INPUT OF A $',A8,
     *          ' GROUP')
 9020 FORMAT(/1X,'**** ERROR, THIS RUN DOES NOT CORRECTLY SPECIFY',
     *          ' THE FULL CI SPACE'/
     *     1X,'CHECK $',A8,' INPUT: NCORE=',I4,' NACT=',I4,' NELS=',I4)
 9030 FORMAT(/1X,'**** ERROR, $ALDET INPUT NCORE=',I4,' NELS=',I4,
     *          ' SZ=',F6.3/
     *       1X,' IS INCONSISTENT WITH $CONTRL INPUT ICH=',I4,
     *          ' MULT=',I4)
 9035 FORMAT(/1X,'***** ERROR, REQUESTED NUMBER OF CI ROOTS=',I5/
     *        1X,'EXCEEDS THE DIMENSION LIMIT FOR NUMBER OF STATES',I5)
 9036 FORMAT(/1X,'**** ERROR, YOUR STATE SELECTED FOR PROPERTIES=',I5/
     *        1X,'EXCEEDS THE NUMBER OF ROOTS YOU REQUESTED=',I5)
 9037 FORMAT(/1X,'**** ERROR, YOUR STATE SELECTED FOR PROPERTIES=',I5/
     *        1X,'EXCEEDS THE DIMENSION LIMIT FOR NUMBER OF STATES',I5)
 9040 FORMAT(/1X,'**** ERROR, WEIGHTS ASSIGNED TO STATES HIGHER',
     *          ' THAN NSTATE=',I5)
 9050 FORMAT(/1X,'**** ERROR, NEGATIVE VALUE FOR -WSTATE- ???')
 9055 FORMAT(1X,'YOUR INPUT ASKS FOR IROOT=',I4,' BUT ONLY',I4,
     *          ' STATE WEIGHTS GIVEN IN WSTATE')
 9061 FORMAT(/1X,'**** WARNING ****'/
     *        1X,'RUNTYP=',A8,' INVOLVING A NUCLEAR GRADIENT'/
     *        1X,'HAS BEEN REQUESTED FOR A RUN AVERAGING OVER',I4,
     *           ' STATES.'/
     *        1X,'THIS MAKES SENSE ONLY IF THE STATES ARE DEGENERATE.'/
     *        1X,'THE RUN IS ALLOWED TO PROCEED BECAUSE USER INPUT'/
     *        1X,'REQUESTS THIS.  PLEASE CHECK THAT THE STATES ARE IN'/
     *        1X,'FACT DEGENERATE, AFTER THEY ARE COMPUTED BELOW.'//)
 9062 FORMAT(/1X,'**** ERROR, RUNTYP=',A8,' REQUIRES NUCLEAR GRADIENT.'/
     *        1X,'THIS IS IMPOSSIBLE WHILE STATE AVERAGING. NAVG=',I5/
     *        1X,'A POSSIBLE RECOURSE, FOR VERY SMALL MOLECULES, IS',
     *           ' TO SELECT'/
     *        1X,'A SPECIFIC STATE ENERGY TO BE NUMERICALLY',
     *           ' DIFFERENTIATED:'/
     *        1X,' $DET    IROOT=NNN'/
     *        1X,' $CONTRL NUMGRD=.TRUE.'//)
 9063 FORMAT(/1X,'**** ERROR, RUNTYP=',A8,' REQUIRES NUCLEAR HESSIAN.'/
     *        1X,'THIS IS IMPOSSIBLE WHILE STATE AVERAGING. NAVG=',I5/
     *        1X,'A POSSIBLE RECOURSE, FOR VERY SMALL MOLECULES, IS',
     *           ' TO SELECT'/
     *        1X,'A SPECIFIC STATE ENERGY TO BE NUMERICALLY',
     *           ' DIFFERENTIATED:'/
     *        1X,' $DET  IROOT=NNN'/
     *        1X,' $HESS METHOD=FULLNUM'//)
 9064 FORMAT(/1X,'**** ERROR, YOUR SA-MCSCF REQUESTS ENERGY',
     *           ' DERIVATIVES, BUT YOU DID'/
     *        1X,'NOT SPECIFY BY -IROOT- WHICH SPECIFIC STATE YOU',
     *           ' WANT.')
C
 9100 FORMAT(/1X,'THE POINT GROUP                  =',3X,A8/
     *       1X,'THE STATE SYMMETRY               =',3X,A8/
     *       1X,'NUMBER OF CORE ORBITALS          =',I5/
     *       1X,'NUMBER OF ACTIVE ORBITALS        =',I5/
     *       1X,'NUMBER OF ALPHA ELECTRONS        =',I5,
     *          ' (',I4,' ACTIVE)'/
     *       1X,'NUMBER OF BETA ELECTRONS         =',I5,
     *          ' (',I4,' ACTIVE)'/
     *       1X,'NUMBER OF OCCUPIED ORBITALS      =',I5)
C
 9160 FORMAT(/1X,'SYMMETRIES FOR THE',I4,' CORE,',I4,' ACTIVE ARE')
      END
