C 25 JUN 07 - TN  - IMPLEMENT AN AUTOMATIC CEEIS METHOD
C
C*MODULE CEEIS   *DECK CEEISCHK
      SUBROUTINE CEEISCHK
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL CEECHK
C
      COMMON /CEEIS0/ ICEEIS
      COMMON /CEEIS3/ NEEDCI, CEECHK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      CEECHK = .FALSE.
      IEOF = 0
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $CEEIS ',IEOF)
      IF (IEOF.NE.0) THEN
        ICEEIS = 0
      ELSE
        ICEEIS = 1
      END IF
      RETURN
      END
C*MODULE CEEIS   *DECK CEEISIN
      SUBROUTINE CEEISIN
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SMRCI,RESTRT,CEMEM,CEMXPN
      LOGICAL FDIRCT,QCORR
C
      PARAMETER (MXATM=2000)
      PARAMETER (MXRT=100)
      PARAMETER (ZERO=0.0D+00)
      PARAMETER (MMX=500,IEXMX=50)
C
CZ    CHARACTER*4 ATM
CZ    CHARACTER*6 BALTYP
CZ    CHARACTER*10 UATOM
CZ    DOUBLE PRECISION MOROKM,MAKEFP,MD,LCCD,NMR,MDTYP,METHMC
CZ    INTEGER DDI_NP,DDI_ME,FLO,CEIL,DIFF
C
CZ    DIMENSION COM(3,MXATM),VMOI(3)
      DIMENSION IEMAT(MMX,IEXMX)
C
CZ    COMMON /INFO  / CDUM(MXATM,3),IZAN(MXATM),NATOMS,IUNTRD,ATM(MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NEE,NAA,NBB,
     *                AN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
CZ    COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST1,JST1,KST1,LST1
CZ
      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 /CEEIS0/ ICEEIS
      COMMON /CEEIS1/ MMIN0,MMIN,IDELTM,SMRCI,M1M2EX(MMX),
     *                XTRTOL(IEXMX),ISCHME,ISTPEX,RESTRT,CEMEM,CEMXPN,
     *                NSEXT(51),IEXPEMT
      COMMON /CEEIS2/ ECONT(IEXMX),EMAT(MMX,IEXMX),ENREF,
     *                MEMCE(MMX,IEXMX),MXPNCE(MMX,IEXMX),
     *                ISTEX(IEXMX),MEX(IEXMX),MMC
CZ    SINCE THIS BLOCK IS ADDED, CEEISIN MUST BE CALLED AFTER FCINPUT
      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
C
C      ----- SET UP NAMELIST $CEEIS TABLES -----
C
      PARAMETER (NNAM=13)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA CEEIS /8HCEEIS   /
      DATA QNAM/8HMMIN    ,8HIDELTM  ,8HSMRCI   ,8HM1M2EX  ,
     *          8HXTRTOL  ,8HISCHME  ,8HISTPEX  ,8HRESTRT  ,
     *          8HCEMEM   ,8HCEMXPN  ,8HENREF   ,8HNSEXT   ,
     *          8HIEXPND  /
      DATA KQNAM /1,1,0,5001,  503,1,1,0,  0,0,3,511, 1/
C
C     ===== CORRELATION ENERGY EXTRAPOLATION BY INTRINSIC SCALING =====
C      PROGRAMMED BY TAKESHI NAGATA, IOWA STATE UNIVERSITY, JUNE 2007.
C     METHOD DEVELOPED BY KLAUS RUEDENBERG, LAIMUTIS BYTAUTAS 2004-2007.
C
C     ----- INITIALIZE VARIABLES FOR NAMELIST $CEEIS -----
C
CZ
      WRITE (6,9000)
C
      NANB = NA + NB
      MREF  = MSTA(NSPACE) - MSTA(1)
      MMIN0 = MAX(MREF,NANB)
      MMIN0 = MMIN0 + NCOR + 1
      MMIN  = MMIN0
      IDELTM  = 1
      ISCHME  = 1
      ISTPEX  = 8
      IEXPEMT = 0
      ENREF  = ZERO ! REFERENCE ENERGY
      SMRCI  = .FALSE.
      RESTRT = .FALSE.
      CEMEM  = .FALSE.
      CEMXPN = .FALSE.
      DO I = 1, MMX
        M1M2EX(I) = 0
      END DO
      DO I = 4, ISTPEX
        XTRTOL(I-3) = 2.0D-04
      END DO
      DO I = ISTPEX+1, IEXMX+3
        XTRTOL(I-3) = ZERO
      END DO
      DO I = 1, 51
        NSEXT(I) = 1
      END DO
CZ
C
C        ----- READ NAMELIST $CEEIS  -----
C
      JRET = 0
      CALL NAMEIO(IR,JRET,CEEIS,NNAM,QNAM,KQNAM,
     *            MMIN,IDELTM,SMRCI,M1M2EX,XTRTOL,
     *            ISCHME,ISTPEX,RESTRT,CEMEM,CEMXPN,
     *            ENREF,NSEXT,IEXPEMT,
     *                       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,0,0,   0,0,0,0,0)
      IF(JRET .EQ. 2) THEN
         IF (MASWRK) WRITE (IW,9005)
         CALL ABRT
      END IF
C
C
C     ---- IF $CEEIS NOT FOUND, RETURN ----
C
      IF(JRET.EQ.1) THEN
         ICEEIS = 0
         RETURN
      END IF
C
      NERR = 0
      IF (MMIN.LT.MMIN0) THEN
         IF (MASWRK) WRITE(IW,9010) 'MMIN', MMIN
         NERR = NERR+1
      END IF
C
      IF (IDELTM.LE.0) THEN
         IF (MASWRK) WRITE(IW,9010) 'IDELTM', IDELTM
         NERR = NERR+1
      END IF
C
      IF (ISTPEX.NE.4.AND.IDELTM.NE.2.AND.SMRCI) THEN
         IF (MASWRK) WRITE(IW,*)
     *   'IF SMRCI IS TRUE, IDELTM AND ISTPEX MUST BE 2 AND 4'
         NERR = NERR+1
      END IF
C
      IF (ISCHME.NE.1.AND.ISCHME.NE.2) THEN
         IF (MASWRK) WRITE(IW,*) 'ISCHME MUST BE 1 OR 2',ISCHME
         NERR = NERR+1
      END IF
C
      IF (ISTPEX.GT.IEXMX) THEN
         IF (MASWRK) WRITE(IW,*)
     *               'ISTPEX MUST BE LESS THAN ',IEXMX
         NERR = NERR+1
      END IF
      IF (ISTPEX.LT.4) THEN
         IF (MASWRK) WRITE(IW,*)
     *               'ISTPEX IS AT LEAST 4, QUADRUPLE'
         NERR = NERR+1
      END IF
C
      IF (ENREF.GT.ZERO) THEN
         IF (MASWRK) WRITE(IW,9015) 'ENREF', ENREF
         NERR = NERR+1
      END IF
C
      DO I = 4, ISTPEX
        IF (XTRTOL(I-3).EQ.ZERO) THEN
          IF (MASWRK) WRITE(IW,9015) 'XTRTOL', XTRTOL(I-3)
          NERR = NERR+1
          GOTO 100
        END IF
      END DO
  100 CONTINUE
C
      IF (CEMEM.OR.CEMXPN) THEN
         IF (MASWRK) WRITE(IW,9015)
     *      'CURRENTLY CEMEM AND CEMXPN ARE NOT AVAILABLE'
         NERR = NERR+1
      END IF
C
CZ
CZ    IEXC: INDICATES WHICH EXCITATION
CZ    IELE: THE NUMBER OF ELEMENTS IN THE EXCITATION
CZ
      IEXC  = 0
      IELE  = 0
      IPREV = 0
      ISTEX(1) = 2
      DO I = 1, 500
C
        M1M2 = M1M2EX(I)
        IF (IELE.EQ.0) THEN
          IEXC = IEXC + 1
          IF (M1M2.EQ.0) THEN
            IELE = 1
          ELSE
            IF (MASWRK) WRITE(IW,*) 'DELIMITER MUST BE 0 ',
     *                              'IN',IEXC+2,'-PULE EXC'
            NERR = NERR+1
          END IF
        ELSEIF (IELE.EQ.1) THEN
          IF (M1M2.EQ.0.OR.M1M2.GE.MMIN) THEN
            IELE  = IELE + 1
            IPREV = M1M2EX(I)
          ELSE
            IF (MASWRK) WRITE(IW,*) M1M2
            IF (MASWRK) WRITE(IW,*) 'FIRST ELEMENT MUST BE MORE THAN ',
     *                         MMIN-1,'OR ZERO IN',IEXC+1,'-PULE EXC'
            NERR = NERR+1
          END IF
        ELSEIF (IELE.EQ.2) THEN
          IF (M1M2.EQ.0) THEN
            IF (IPREV.EQ.0) THEN
              MEX(IEXC) = IELE
              ISTEX(IEXC+1)  = I + 2
              IELE  = 0
            ELSE
              IF (MASWRK) WRITE(IW,*) 'VALUE IN ELEMENT ',I,
     *           ' IS 0 THOUGH PREV VALUE IS OTHER THAN 0'
              NERR = NERR+1
            ENDIF
          ELSEIF (M1M2.GE.MMIN) THEN
            ! FOR IELE=2, M1M2EX > MMIN
            IF (IPREV.EQ.0) THEN
              IF (MASWRK) WRITE(IW,*) 'VALUE IN ELEMENT ',I,
     *           ' IS OTHER THAN 0 THOUGH PREV VALUE IS 0'
              NERR = NERR+1
            ELSEIF (M1M2.GT.IPREV) THEN
              IELE  = IELE + 1
              IPREV = M1M2
            ELSE
              IF (MASWRK) WRITE(IW,*) 'ELEMENTS MUST BE ',
     *           'INCREMENTAL NUMBERS'
              NERR = NERR+1
            ENDIF
          ELSEIF (M1M2.LT.0) THEN
            IF (IPREV.EQ.0) THEN
              IF (MASWRK) WRITE(IW,*) 'VALUE IN ELEMENT ',I,
     *           ' IS NEGATIVE THOUGH PREV VALUE IS 0'
              NERR = NERR+1
            ELSEIF (-M1M2.GT.IPREV) THEN
              IELE  = IELE + 1
              IPREV = M1M2
            ELSE
              IF (MASWRK) WRITE(IW,*) 'ABSOLUTE VALUE IN ELEMENT ',
     *           I, ' IS SMALLER THAN THE PREVIOUS VALUE'
              NERR = NERR+1
            ENDIF
          ENDIF
        ELSEIF (IELE.GE.3) THEN
          IF (M1M2.EQ.0) THEN
            MEX(IEXC) = IELE - 1
            IEXC = IEXC + 1
            ISTEX(IEXC) = I + 1
            IELE  = 1
            IPREV = 0
          ELSEIF (M1M2.GT.MMIN) THEN
            IELE  = IELE + 1
            IPREV = M1M2
          ELSEIF (M1M2.LT.0) THEN
            IF (IPREV.GT.0.AND.-M1M2.GT.IPREV) THEN
              IELE  = IELE + 1
              IPREV = M1M2
            ELSE
              IF (MASWRK) WRITE(IW,*) 'VALUE IN ELEMENT ',I-1,
     *           ' MUST BE CHANGED'
              NERR = NERR+1
            ENDIF
          ELSE
            IF (MASWRK) WRITE(IW,*) 'VALUE IN ELEMENT ',I,
     *         ' IS LESS THAN ', MMIN
            NERR = NERR+1
          ENDIF
        ENDIF
C
        IF (IEXC+2.GT.ISTPEX) GOTO 101
      END DO
  101 CONTINUE
C
CZ
CZ    PUT RISTRICTION FOR RECURSIVE CEEIS
CZ    EX: M1M2EX(1)=0,0,0,0,0,0,0,0,0,0 --> OK
CZ        M1M2EX(1)=0,?,?,0,0,0,0,0,0,0 --> OK
CZ        M1M2EX(1)=0,?,?,0,0,0,0,%,%,0 --> KILLED
CZ        M1M2EX(1)=0,0,0,0,%,%,0,0,0,0 --> KILLED
CZ        ? AND % ARE NUMBERS
CZ
      ICUNT = 0
      NCUNT = 0
      DO III = 2, ISTPEX - 2
        IF (III.EQ.2.OR.ISCHME.EQ.1) THEN
          NCUNT = NCUNT + 1
        ELSEIF (III.GE.3) THEN
          IF (MOD(III,2).EQ.0) THEN
            NCUNT = NCUNT + 1
          ELSE
            GOTO 105
          END IF
        END IF
        IF (M1M2EX(ISTEX(III)).EQ.0) THEN
          ICUNT = ICUNT + 1
        ENDIF
  105   CONTINUE
      END DO
      IF (ICUNT.GT.0.AND.ICUNT.NE.NCUNT) THEN
        WRITE(6,*) 'IF YOU USE RECURSIVE CEEIS, M1M2EX(1) ',
     *             'MUST BE ALL ZEROS EXCEPT FOR TRIPLE EXC'
        WRITE(6,*) 'ISTPEX MUST BE CHOSEN EQUAL TO HIGHEST ',
     *             'SPECIFIED EXCITATION LEVEL IN M1M2EX(1)'
        NERR = NERR + 1
      ENDIF
      IF (ICUNT.EQ.0.AND.SMRCI) THEN
        WRITE(6,*) 'SMRCI MUST BE RECURSIVE CEEIS'
        NERR = NERR + 1
      ENDIF
      DO I = 1, 51
        IF (NSEXT(I).NE.1.AND.NSEXT(I).NE.0) THEN
          WRITE(6,*) 'NSEXT(1) MUST BE 0 OR 1'
          NERR = NERR + 1
        END IF
      END DO
      IF (IEXPEMT.LT.0) THEN
        WRITE(6,*) 'IEXPEMT MUST BE POSITIVE'
        NERR = NERR + 1
      END IF
CZ
CZ    READ $CEEIS GROUP
CZ
CZ    THE NUMBER OF CORRELATING ORBITALS
CZ    MMC = NQMT - NCOR - (NA+NB)/2
      MMC = NQMT
      IF (RESTRT) CALL RDCEEIS
CZ
      DO J = 1, ISTPEX - 1
        DO I = 1, MMC - MMIN + 1
          IEMAT(I,J) = 0
          IF (EMAT(I,J).LT.ZERO) THEN
            IEMAT(I,J) = 1
          END IF
        END DO
      END DO
CZ
CZ
CZ    I-PLE EXC MUST CONTAINE SOME OF ORBITALS FOR I+1(I+2)-PLE EXC
CZ
      DO III = 1, ISTPEX-3
        IEX1 = ISTEX(III)
        IF (III.GE.2.AND.ISCHME.EQ.2) THEN
          IF (MOD(III,2).EQ.0) THEN
            NEX  = III + 2
            IEX2 = ISTEX(III+2)
          ELSE
            GOTO 104
          END IF
        ELSE
          NEX  = III + 1
          IEX2 = ISTEX(III+1)
        END IF
        IF (M1M2EX(IEX1).NE.0.AND.M1M2EX(IEX2).NE.0) THEN
          I = 0
  102     CONTINUE
            I = I + 1
            IF (M1M2EX(IEX2+I).LT.0) THEN
              ISTART =  M1M2EX(IEX2+I-1)
              IEND   = -M1M2EX(IEX2+I)
              I = I + 1
            ELSE
              ISTART =  M1M2EX(IEX2+I-1)
              IEND   =  ISTART
            ENDIF
C
            DO J = ISTART, IEND
              IOK = 0
CZ RESTART CHECK
              DO KJ = 1, MMC - MMIN + 1
                IF (J.EQ.KJ+MMIN-1) THEN
                  IF (IEMAT(KJ,III+1).EQ.1) GOTO 106
                END IF
              END DO
CZ
              KK  = 0
  103         CONTINUE
                KK = KK + 1
                IF (M1M2EX(IEX1+KK).LT.0) THEN
                  KSTART =  M1M2EX(IEX1+KK-1)
                  KEND   = -M1M2EX(IEX1+KK)
                  KK = KK + 1
                ELSE
                  KSTART =  M1M2EX(IEX1+KK-1)
                  KEND   =  KSTART
                ENDIF
                DO L = KSTART, KEND
                  IF (L.EQ.J) GOTO 106
                END DO
              IF (KK.LT.MEX(III)) GOTO 103
              IOK = IOK + 1
  106         CONTINUE
              IF (IOK.NE.0) THEN
                WRITE (6,*) 'FOR EXC ',III+3,' OR ',III+4,
     *                      'REDUNDANT NUMBER ',J,
     *                      ' IS INCLUDED IN M1M2EX(1)'
                NERR = NERR + 1
              END IF
            END DO
          IF (I.LT.MEX(NEX)) GOTO 102
        ENDIF
  104   CONTINUE
      END DO
CZ
C
C
C     ----- ONE OR MORE ERRORS DETECTED, BLOW THE JOB AWAY -----
C
      IF(NERR.GT.0) THEN
         IF (MASWRK) WRITE(IW,9200) NERR
         CALL ABRT
      END IF
CZ
C
C
C        PRINT $CEEIS
C
      IF (MASWRK) THEN
        WRITE (IW,9520) MMIN  , IDELTM, SMRCI ,  ISCHME,
     *                  ISTPEX, RESTRT, CEMEM ,  CEMXPN,
     *                  ENREF,  IEXPEMT
        WRITE (IW,9525)
        DO I = 1, ISTPEX-2
          WRITE (IW,9530) I+2,(M1M2EX(ISTEX(I)+J),J=0,MEX(I)-1)
        END DO
        WRITE (IW,9540) (XTRTOL(I-3),I=4,ISTPEX)
        WRITE (IW,9600)
        IF (NCOR.NE.0) THEN
          WRITE (IW,9650) 1, NCOR, 2*NCOR
        END IF
        DO I = 1, NSPACE-1
          IF (NSEXT(I).EQ.1) THEN
            WRITE (IW,9700) MSTA(I),MSTA(I+1)-1,MINI(I),MAXI(I)
          ELSE
            WRITE (IW,9750) MSTA(I),MSTA(I+1)-1,MINI(I),MAXI(I)
          END IF
        END DO
        WRITE (IW,9800) MSTA(NSPACE),MMC,MINI(I),MAXI(I)
        WRITE (IW,9850) NA+NB+NCOR*2
        WRITE (IW,9550)
      END IF
C
 9000 FORMAT(/3X,57(1H-)/
     *  3X,'CEEIS CALCULATION'/
     *  3X,'TAKESHI NAGATA, LAIMUTIS BYTAUTAS, KLAUS RUEDENBERG'/
     *  3X,'AMES LABORATORY USDOE, IOWA STATE UNIVERSITY, 2007'/
     *  3X,'REFERENCES:  L. BYTAUTAS AND K.RUEDENBERG, J. CHEM. PHYS.'/
     *  3X,'121, 10905 (2004); 121, 10919 (2004); 122, 154110 (2005).'/
     *  3X, 57(1H-))
 9005 FORMAT(1X,'TYPING ERROR IN $CONTRL INPUT - CHECK NEAR $ MARKER')
 9010 FORMAT(/1X,'ERROR: $CEEIS KEYWORD ',A6,
     *          ' WAS GIVEN AN ILLEGAL VALUE ',I8,'.'/)
 9015 FORMAT(/1X,'ERROR: $CONTRL KEYWORD ',A6,
     *          ' WAS GIVEN AN ILLEGAL VALUE ',F10.8,'.'/)
 9200 FORMAT(/1X,'             *** ERROR(S) DETECTED ***'/
     *        1X,'YOUR $CEEIS  INPUT HAS AT LEAST ONE SPELLING OR',
     *           ' LOGIC MISTAKE.',I4/
     *        1X,'PLEASE REVIEW THE REASON(S) JUST LISTED, AND TRY',
     *           ' YOUR RUN AGAIN.'/)
 9520 FORMAT(/5X,'$CEEIS  OPTIONS'/5X,15(1H-)/
     * 4X,'MMIN  =',I8,3X,'IDELTM=',I8,3X,'SMRCI =',L8,3X,'ISCHME=',I8/
     * 4X,'ISTPEX=',I8,3X,'RESTRT=',L8,3X,'CEMEM =',L8,3X,'CEMXPN=',L8/
     * 4X,'ENREF =',F16.8,3X,'IEXPEMT=',I8)
 9525 FORMAT(/4X,50(1H-)/
     *       4X,'ENERGIES REQUESTED FOR'/
     *       4X,'EXCITATION LEVELS     HIGHEST ORBITALS USED'/
     *       4X,50(1H-))
 9530 FORMAT(11X,'X=',I2,20X,'M=',8I5)
 9540 FORMAT(/5X,'THRESHOLDS FOR EXCITATION ENERGIES'/
     *        5X, 8(1P,E8.1))
C9545 FORMAT(/5X,'ORMAS SPACE ACTIVATION',8I5)
 9550 FORMAT(3X,70(1H=))
 9600 FORMAT(/4X,'ORMAS SEGMENTATION PATTERN OF REFERENCE FUNCTION')
 9650 FORMAT(4X,'     CORE ORBITALS',I3,'-',I3,2X,' OCCUPATIONS',I3)
 9700 FORMAT(4X,'REFERENCE ORBITALS',I3,'-',I3,2X,' OCCUPATIONS',
     *          I3,'-',I3,2X,'EXCITATION-ACTIVE')
 9750 FORMAT(4X,'REFERENCE ORBITALS',I3,'-',I3,2X,' OCCUPATIONS',
     *          I3,'-',I3,2X,'EXCITATION-INACTIVE')
 9800 FORMAT(4X,'  VIRTUAL ORBITALS',I3,'-',I3,2X,' OCCUPATIONS',
     *          I3,'-',I3)
 9850 FORMAT(4X,'TOTAL NUMBER OF ELECTRONS',I4)
      END
C*MODULE CEEIS   *DECK RDCEEIS
      SUBROUTINE RDCEEIS
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
CZ    PARAMETER (MXAO=8192)
C
CZ    LOGICAL DBG,MIX,GOPARR,DSKWRK,MASWRK,MFRZ,MOIDON,EDCOMP,
CZ   *        DIPDCM,DEPRNT,QADDCM,ZDO,POLDCM,POLANG,POLAPP,
CZ   *        KMIDPT,POLDYN
      LOGICAL GOPARR,DSKWRK,MASWRK,SMRCI,RESTRT,CEMEM,CEMXPN
      CHARACTER*8 STCEDATA
CZ    DOUBLE PRECISION MOROKM
C
CZ    DIMENSION IWRK(L1),EIG(L1),SCR(L1,8),S(L2),Q(L3),VEC(L1,L1)
C
CZ    PARAMETER (MXATM=2000)
CZ    PARAMETER (MXRT=100)
      PARAMETER (MMX=500,IEXMX=50)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
CZ    COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NEE,NAA,NBB,
CZ   *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
CZ    COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
CZ   *                NCO,NSETO,NOPEN,NPAIR,NMOGVB,NCONF(MXAO),NHAM
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
CZ    COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
CZ
CZ    I NEED NACT AND NA
CZ    COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,S,SZ,
CZ   *                GRPDET,STSYM,GLIST,
CZ   *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORB,
CZ   *                NA,NB,K,KST,IROOT,IPURES,MAXW1,NITER,MAXP,NCI,
CZ   *                IGPDET,KSTSYM,NFTGCI
      COMMON /CEEIS1/ MMIN0,MMIN,IDELTM,SMRCI,M1M2EX(MMX),
     *                XTRTOL(IEXMX),ISCHME,ISTPEX,RESTRT,CEMEM,CEMXPN,
     *                NSEXT(51),IEXPEMT
      COMMON /CEEIS2/ ECONT(IEXMX),EMAT(MMX,IEXMX),ENREF,
     *                MEMCE(MMX,IEXMX),MXPNCE(MMX,IEXMX),
     *                ISTEX(IEXMX),MEX(IEXMX),MMC
C
C     THIS ROUTINE MUST BE CALLED AFTER READING $CIDET GROUP
C
      IZERO = 0
      IONE  = 1
      ITWO  = 2
C
      STCEDATA = ' $CEDATA'
C
      IF(MMC.EQ.0.AND.MASWRK) WRITE(IW,9080)
      IF(MMC.EQ.0) CALL ABRT
C
C     POSITION THE INPUT FILE TO $CEDATA
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,STCEDATA,IEOF)
      IF (IEOF.EQ.0) GO TO 200
C
C               $CEDATA CARD NOT FOUND, KILL JOB
C
      IF (MASWRK) WRITE(IW,9020)
      CALL ABRT
      STOP
C
C               $CEDATA CARD FOUND, READ THE ORBITALS
C
C
  200 CONTINUE
CZ
      IF (MASWRK) THEN
CZ    IC = 0
      IC  = MMIN - 1
      III = MMC - MMIN + 1
      JJJ = ISTPEX - IEXPEMT - 1
      DO I = 1, III
         IMAX = 0
         IC = IC + 1
  240    CONTINUE
            IMIN = IMAX+1
            IMAX = IMAX+5
            IF(IMAX .GT. JJJ) IMAX = JJJ
            READ(IR,9040,END=300,ERR=300) ICC, (EMAT(I,J),J=IMIN,IMAX)
            IF(ICC.EQ.IC) GO TO 260
               WRITE(IW,9060) IC, ICC
               IF (GOPARR) CALL DDI_BCAST(350,'I',ITWO,1,MASTER)
               CALL ABRT
               STOP
  260       CONTINUE
         IF(IMAX .LT. JJJ) GO TO 240
      END DO
CZ
C
      IF (GOPARR) CALL DDI_BCAST(350,'I',IZERO,1,MASTER)
C
C SLAVE WORK
C
      ELSE
         IF (GOPARR) CALL DDI_BCAST(350,'I',IEND,1,MASTER)
         IF (IEND.EQ.IONE) GO TO 300
         IF (IEND.EQ.ITWO) CALL ABRT
      END IF
C
C  GIVE VECTORS TO ALL PROCESSES
C
CZ    IF (GOPARR) CALL DDI_BCAST(351,'F',VEC,NMOS*L1,MASTER)
      IF (GOPARR) CALL DDI_BCAST(351,'F',EMAT,III*JJJ,MASTER)
CZ
      RETURN
CZ
C
C        PREMATURE END OF $CEDATA INPUT ENCOUNTERED
C
  300 CONTINUE
      IF (MASWRK) THEN
         IF (GOPARR) CALL DDI_BCAST(350,'I',IONE,1,MASTER)
         WRITE(IW,9100) ICC
      END IF
      CALL ABRT
 9020 FORMAT(32H NO $CEDATA INPUT FOUND - STOP  )
 9040 FORMAT(I2,5F16.8)
 9060 FORMAT(' *** ERROR READING THE ORMAS ENERGY MATRIX AT',I4/
     *       '     THE INPUT VALUES WERE',I4)
 9080 FORMAT(1X,'NACT AND NELS ARE WRONG IN $CIDET')
 9100 FORMAT(' *** ERROR: PREMATURE END OF CEDATA INPUT ENCOUNTERED'/
     *       '            LOOKING FOR ELEMENT',I4)
      END
C*MODULE CEEIS    *DECK CEEIS
C     ------------------------------
      SUBROUTINE CEEIS(NRNFG,NPFLG)
      ! THIS ROUTINE IS TAKEN FROM ORDET IN ORMAS1.SRC
C     ------------------------------
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION NRNFG(10),NPFLG(10)
CZ
C
      LOGICAL GOPARR,DSKWRK,MASWRK,DOEXCH,
     *        DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,DOVVVO,DOVVVV,DOCORE
      LOGICAL FDIRCT,QCORR
      LOGICAL SMRCI,RESTRT,CEMEM,CEMXPN
      LOGICAL CEECHK
C
CZ    PARAMETER (MXATM=2000, NMO=500)
      PARAMETER (NMO=500)
      PARAMETER (MXRT=100)
      PARAMETER (ZERO=0.0D+00)
      PARAMETER (MMX=500,IEXMX=50)
CZ
      CHARACTER*14 CDELE(3), CDELE1(8), CDELE2(8)
      DIMENSION COEF(3), XX(NMO), YY(NMO), ZZ(NMO)
      DIMENSION DEMAT(MMX,IEXMX)
      DIMENSION MAXISV(51), MINISV(51)
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
C
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,STOT,SSQUAR,ECORE,ESCF,EERD,
     *                E1,E2,VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
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
C
      COMMON /IOFILE/ IR,IW,IP,IS,IJKT,IDAF,NAV,IODA(950)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /CEEIS1/ MMIN0,MMIN,IDELTM,SMRCI,M1M2EX(MMX),
     *                XTRTOL(IEXMX),ISCHME,ISTPEX,RESTRT,CEMEM,CEMXPN,
     *                NSEXT(51),IEXPEMT
      COMMON /CEEIS2/ ECONT(IEXMX),EMAT(MMX,IEXMX),ENREF,
     *                MEMCE(MMX,IEXMX),MXPNCE(MMX,IEXMX),
     *                ISTEX(IEXMX),MEX(IEXMX),MMC
      COMMON /CEEIS3/ NEEDCI, CEECHK
      COMMON /FMPARM/ LTOP,LOFFS,LENHI,LOCMEM,MEMLIM,MEMOK
      COMMON /CIFILS/ NFT11,NFT12,NFT13,NFT14,NFT15,NFT16,IDAF20,NEMEMX
C
      DATA CDELE /'    DEMAT(M,2)','    DEMAT(M,3)','    DEMAT(M,4)'/
      DATA CDELE1/'    DEMAT(M,5)','    DEMAT(M,6)','    DEMAT(M,7)',
     *            '    DEMAT(M,8)','    DEMAT(M,9)','   DEMAT(M,10)',
     *            '   DEMAT(M,11)','   DEMAT(M,12)'/
      DATA CDELE2/'              ','   DEMAT(M,56)','              ',
     *            '   DEMAT(M,78)','              ','  DEMAT(M,910)',
     *            '              ',' DEMAT(M,1112)'/
CZ
C
      DATA CIDET/8HCIDET    /
      DATA ORMAS/8HORMAS    /
      DATA RNONE/8HNONE     /
CZ
CZ    SCHEME  I: YY = COEF(1) + COEF(2)*XX
CZ    SCHEME II: ZZ = COEF(1)*XX + COEF(2)*YY + COEF(3)
CZ
C
C        DRIVER FOR FULL CLASS CI CALCULATIONS...
C
      DO III = 1, IEXMX
        ECONT(III) = ZERO
        DO JJJ = 1, MMX
           EMAT(JJJ,III) = ZERO
          DEMAT(JJJ,III) = ZERO
        END DO
      END DO
      DO III = 1, NMO
        XX(III) = ZERO
        YY(III) = ZERO
        ZZ(III) = ZERO
      END DO
C
C        ----- READ INPUT DEFINING THE FULL CLASS CI DIMENSIONS -----
C
      CALL FCINPUT(NPFLG(1),CIDET,ORMAS)
      CALL CEEISIN
CZ
      DO KK = 1, NSPACE
        MINISV(KK) = MINI(KK)
        MAXISV(KK) = MAXI(KK)
      END DO
CZ
CZ    IF THE REF ENERGY IS NOT GIVEN, DO THE CALCULATION HERE
CZ
      IF (ENREF.EQ.ZERO) THEN
        IF (MASWRK) THEN
          WRITE(IW,'(/3X,40(1H-))')
          WRITE(IW,'(3X,A22,A16)') 'ORMAS CALCULATION FOR ',
     *                           'REFERENCE ENERGY'
        END IF
        CALL ORDET(NRNFG,NPFLG)
        CALL SEQCLO(IJKT,'DELETE')
        CALL SEQCLO(NFT12,'DELETE')
        ENREF = ETOT
        IF (MASWRK) THEN
          WRITE(IW,'(3X,A6,F16.10)') 'ENREF=', ENREF
          WRITE(IW,'(3X,40(1H-)/)')
        END IF
      END IF
CZ
      NANB = (NA + NB)/2
      NORTOT = NANB + NCORSV
C
      SIGMAY = 1.0D+00
CZ
CZ    OUTER LOOP: EXCITATIONS
CZ
      DO III = 2, ISTPEX
C
        IF (ISCHME.EQ.2.AND.III.GE.5.AND.MOD(III,2).NE.0)
     *  GOTO 4000
CZ
        IF (MASWRK) THEN
          WRITE(IW,'(/3X,44(1H=))')
          WRITE(IW,'(3X,A33,I4)')
     *            'CALCULATIONS FOR EXCITATION LEVEL', III
          WRITE(IW,'(3X,44(1H=)/)')
        ENDIF
CZ
CZ
        IEXTRA = 0
        MINTMP = MMIN
        IF (III.EQ.2.OR.III.EQ.3) THEN
          IF (III.EQ.2) IEXTRA = 1
          IF (M1M2EX(ISTEX(1)).EQ.0) THEN
            IEXTRA = 1
            IF (M1M2EX(ISTEX(2)).EQ.0) THEN
              MAXTMP = MINTMP + 19
              IF (MMC.LE.MAXTMP) MAXTMP = MMC
              IPOINT = 0
              IEND   = 0
            ELSE
              IPOINT = ISTEX(2)
              IEND   = ISTEX(2) + MEX(2) - 1
            END IF
          ELSE
            IPOINT = ISTEX(1)
            IEND   = ISTEX(1) + MEX(1) - 1
          END IF
        ELSEIF (III.GE.4) THEN
          IF (M1M2EX(ISTEX(III-2)).EQ.0) THEN
            IF (IDELTM.EQ.1) THEN
              MAXTMP = MINTMP + 4
            ELSE
              MAXTMP = MINTMP + 9
            ENDIF
            IPOINT = 0
            IEND   = 0
          ELSE
            IPOINT = ISTEX(III-2)
            IEND   = ISTEX(III-2) + MEX(III-2) - 1
          ENDIF
        ENDIF
C
CZ
CZ      RANGE CYCLE OR LOOP CONTROLED BY IPOINT
CZ
        ICUNT = 0
        LCUNT = 0
  100   CONTINUE
C
          IF (IPOINT.EQ.0) THEN
            ICUNT = ICUNT + 1
            IF (III.GE.4.AND.ICUNT.GT.1) THEN
              XX(1) = XX(LCUNT-3)
              XX(2) = XX(LCUNT-2)
              XX(3) = XX(LCUNT-1)
              YY(1) = YY(LCUNT-3)
              YY(2) = YY(LCUNT-2)
              YY(3) = YY(LCUNT-1)
              IF (III.GE.5.AND.ISCHME.EQ.2) THEN
                ZZ(1) = ZZ(LCUNT-3)
                ZZ(2) = ZZ(LCUNT-2)
                ZZ(3) = ZZ(LCUNT-1)
              ENDIF
              LCUNT = 3
            ELSE
              LCUNT = 0
            ENDIF
          ELSEIF (IPOINT.LE.IEND) THEN
            IF (M1M2EX(IPOINT+1).LT.0) THEN
              MINTMP =  M1M2EX(IPOINT)
              MAXTMP = -M1M2EX(IPOINT+1)
              IPOINT =  IPOINT + 2
            ELSE
              MINTMP =  M1M2EX(IPOINT)
              MAXTMP =  MINTMP
              IPOINT =  IPOINT + 1
            ENDIF
          ENDIF
C
CZ
CZ        MM DO LOOP STARTS
CZ
          MINTMP = MINTMP - NORTOT
          MAXTMP = MAXTMP - NORTOT
          DO JJJ = MINTMP, MAXTMP, IDELTM
C
            LLL = JJJ - MMIN + NORTOT + 1
            IF (RESTRT.AND.EMAT(LLL,III-1).LT.ZERO) GOTO 5000
C
            NACT = NANB + JJJ
C
            NORB = NCORSV + NACT
            MSTA(NSPACE+1) = NACT + NCORSV + 1
            MNUM(NSPACE) = MSTA(NSPACE+1)-MSTA(NSPACE)
C
            IF (MASWRK) WRITE(IW,9000) III,MNUM(NSPACE),NCORSV,NACT
C
            DO KKK = 1, NSPACE - 1
              IF (NSEXT(KKK).EQ.1) THEN
                NELSCAP = MNUM(KKK)*2
                MINI(KKK) = MINISV(KKK) - III
                MAXI(KKK) = MAXISV(KKK) + III
                IF (MINI(KKK).LT.0) MINI(KKK) = 0
                IF (MAXI(KKK).GT.NA+NB) MAXI(KKK) = NA + NB
                IF (MAXI(KKK).GT.NELSCAP) MAXI(KKK) = NELSCAP
              ELSE
                MINI(KKK) = MINISV(KKK)
                MAXI(KKK) = MAXISV(KKK)
              END IF
            END DO
            IF (NSEXT(NSPACE).EQ.1) THEN
              MINI(NSPACE) = MINISV(NSPACE) - III
              MAXI(NSPACE) = MAXISV(NSPACE) + III
            ELSE
              MINI(NSPACE) = MINISV(NSPACE)
              MAXI(NSPACE) = MAXISV(NSPACE)
            END IF
            IF (MINI(NSPACE).LT.0) MINI(NSPACE)=0
            NELSCAP = MNUM(NSPACE)*2
            IF (MAXI(NSPACE).GT.NELSCAP) MAXI(NSPACE) = NELSCAP
            WRITE(IW,9050) (MNUM(I),I=1,NSPACE)
            WRITE(IW,9060) (MINI(I),I=1,NSPACE)
            WRITE(IW,9070) (MAXI(I),I=1,NSPACE)
C
C
CZ          MAXI, MNUM ETC MUST BE CHECKED IN EACH CEEIS CYCLE
            CALL FCCHECK(IW,.FALSE.,
     *             NSPACE,MNUM,MINI,MAXI,IAMI,IAMA,IBMI,IBMA,
     *             NA,NB)
CZ
C
C           ----- INTEGRAL TRANSFORMATION -----
C
            DDITRF=GOPARR
            DOOOOO=.TRUE.
            DOVOOO=.FALSE.
            DOVVOO=.FALSE.
            DOVOVO=.FALSE.
            DOVVVO=.FALSE.
            DOVVVV=.FALSE.
            DOCORE=.TRUE.
            DOEXCH=SCFTYP.EQ.RNONE
            CALL TRFMCX(NPFLG(2),NCORSV,NORB,NORB,.FALSE.,DOEXCH,
     *                  DDITRF,DOOOOO,DOVOOO,DOVVOO,DOVOVO,
     *                  DOVVVO,DOVVVV,DOCORE)
C
            CEECHK = .TRUE.
            MXPTMP = MAXP
 6000       CONTINUE
C
C           ----- DIRECT ORMAS CI CALCULATION -----
C
            CALL DEFCCI(NPFLG(3),.FALSE.,DDITRF,
     *                  NRNFG(5),NRNFG(6),NPFLG(5),NPFLG(6))
C
            IF (NEEDCI*1.01D+00.GE.MEMLIM) THEN
              MAXP = MAXP - 1
              IF (MAXP.LT.2) THEN
                WRITE(IW,9220) MAXP, NEEDCI, MEMLIM
                CALL ABRT
              ENDIF
              GOTO 6000
            ELSE
              IF (CEECHK) THEN
                CEECHK = .FALSE.
                WRITE(IW,9150) STSYM,GRPDET,SZ,NCI
                GOTO 6000
              ENDIF
              WRITE(IW,*) 'MAXP=',MAXP
              MAXP = MXPTMP
            ENDIF
C
C
CZ          THIS IS TAKEN FROM SUBROUTINE FCINPUT
            IF (QCORR) THEN
              IF (NSPACE.GT.3.OR.NSPACE.LT.2) GO TO 2000
              IF (MAXI(NSPACE).NE.2) GO TO 2000
              IF (MINI(NSPACE).NE.0) GO TO 2000
              IF (MAXI(1)-MINI(1).NE.2) GO TO 2000
              GO TO 3000
 2000         QCORR=.FALSE.
 3000         CONTINUE
            ENDIF
CZ
            IF (QCORR) THEN
C
              WRITE(IW,9010)
              EMRCISD = ETOT
              WRITE(IW,9020) EMRCISD
              WRITE(IW,9025) C0SQ
C
              NSPACE = NSPACE - 1
              IF (NSPACE.GT.1) THEN
                MINI(NSPACE) = MINI(NSPACE) + 2
                MAXI(NSPACE) = MAXI(NSPACE) - 2
              ENDIF
CZ
              MAXMAX       = MAXI(1)
CZ
              MINI(1)      = MINI(1) + 2
              MAXI(1)      = MINI(1)
C
CZ            WRITE(IW,9030)
CZ            WRITE(IW,9040) NSPACE,(MSTA(I),I=1,NSPACE)
CZ            WRITE(IW,9050) (MNUM(I),I=1,NSPACE)
CZ            WRITE(IW,9060) (MINI(I),I=1,NSPACE)
CZ            WRITE(IW,9070) (MAXI(I),I=1,NSPACE)
C
              CALL FCCHECK(IW,.FALSE.,
     *               NSPACE,MNUM,MINI,MAXI,IAMI,IAMA,IBMI,IBMA,
     *               NA,NB)
C
              CALL DEFCCI(-5,.FALSE.,DDITRF,
     *                    0,0,NPFLG(5),NPFLG(6))
C
CZ            IF FORTH AND FIFTH ARGUES ARE BOTH 0,
CZ            MSTA VALUES ARE NOT RECOVERED!
              DO IJK = 1, NSPACE+1
                MSTA(IJK) = MSTA(IJK) + NCORSV
              END DO
C
C
              EREF0 = ETOT
              WRITE(IW,9080) EREF0
              WRITE(IW,9090)
C
              ELOW = EREF0 - EMRCISD
              EQ   = ELOW*(1.0D+00 - C0SQ)/C0SQ
              ESDQ = EMRCISD - EQ
              WRITE(IW,9100) ESDQ
              WRITE(IW,9110)
              WRITE(IW,9120)
C
              ETOT = EMRCISD
C
CZ
CZ            A TRICKY PART FOR CEEIS
CZ
              IF (NSPACE.GT.1) THEN
                 MINI(NSPACE) = MINI(NSPACE) - 2
                 MAXI(NSPACE) = MAXI(NSPACE) + 2
              ENDIF
              NSPACE = NSPACE + 1
              MINI(1)      = MINI(1) - 2
              MAXI(1)      = MAXMAX
CZ
CZ            END OF THE TRICKY PART
CZ
C
            ENDIF
CZ
            EMAT(LLL,III-1) = ETOT
            DO I1=1,K
              WRITE(IW,9004) I1,ESTATE(I1),SPINS(I1),SZ,STSYM
            END DO
            WRITE (6,9005) NACT+NCORSV, III, ETOT
CZ
CZ  5000: SKIP ORMAS CALCULATION
CZ
 5000       CONTINUE
            LCUNT = LCUNT + 1
            IF (III.EQ.2) THEN
              DEMAT(LLL,1) = EMAT(LLL,1) - ENREF
            ELSEIF (III.EQ.3) THEN
              DEMAT(LLL,2) = EMAT(LLL,2) - EMAT(LLL,1)
              IF (M1M2EX(ISTEX(1)).GT.0) THEN
                XX(LCUNT)  = DEMAT(LLL,1)
                YY(LCUNT)  = DEMAT(LLL,2)
              END IF
            ELSEIF (III.EQ.4) THEN
              DEMAT(LLL,1) = EMAT(LLL,1) - ENREF
              DEMAT(LLL,3) = EMAT(LLL,3) - EMAT(LLL,2)
              XX(LCUNT)    = DEMAT(LLL,1)
              YY(LCUNT)    = DEMAT(LLL,3)
            ELSEIF (III.GE.5) THEN
              IF (ISCHME.EQ.1) THEN
                DEMAT(LLL,III-3) = EMAT(LLL,III-3) - EMAT(LLL,III-4)
                DEMAT(LLL,III-1) = EMAT(LLL,III-1) - EMAT(LLL,III-2)
                XX(LCUNT)        = DEMAT(LLL,III-3)
                YY(LCUNT)        = DEMAT(LLL,III-1)
              ELSE
                DEMAT(LLL,III-1) = EMAT(LLL,III-1) - EMAT(LLL,III-3)
                XX(LCUNT)        = DEMAT(LLL,1)
                YY(LCUNT)        = DEMAT(LLL,2)
                ZZ(LCUNT)        = DEMAT(LLL,III-1)
              ENDIF
            ENDIF
C
            CALL SEQCLO(IJKT,'DELETE')
            CALL SEQCLO(NFT12,'DELETE')
C
CZ
CZ          WRITE OUT $CEDATA
CZ
            IF (MASWRK) THEN
              WRITE(IP,9230) '$CEDATA'
              IC = MMIN - 1
              IJK = MMC - MMIN + 1
              JJ2 = ISTPEX - 1
              DO I1 = 1, IJK
                IMAX = 0
                IC = IC + 1
  230           CONTINUE
                  IMIN = IMAX+1
                  IMAX = IMAX+5
                  IF(IMAX .GT. JJ2) IMAX = JJ2
                  WRITE(IP,9210) IC, (EMAT(I1,J1),J1=IMIN,IMAX)
                IF(IMAX .LT. JJ2) GO TO 230
              END DO
              WRITE(IP,9230) '$END'
            END IF
CZ
          END DO
CZ
CZ        END OF ACTIVE ORBITAL LOOP
CZ
C
        IF (IPOINT.LE.IEND) THEN
          GOTO 100
        ELSEIF (IPOINT.GT.IEND) THEN
          IF ((III.EQ.2.OR.III.EQ.3).AND.IEXTRA.EQ.1) THEN
            IEXTRA = 0
            MINTMP = MMC
            MAXTMP = MMC
            GOTO 100
          END IF
        END IF
C
CZ
CZ      EXTRAPOLATION PART
CZ
        IF (III.EQ.2) THEN
          ECONT(1) = EMAT(LLL,1) - ENREF
          IF (MASWRK) THEN
            WRITE(IW,9160) III
            WRITE(IW,9170) ECONT(1)
            WRITE(IW,9190)
          END IF
        ELSEIF (III.EQ.3.AND.M1M2EX(ISTEX(III-2)).EQ.0) THEN
          ECONT(2) = EMAT(LLL,2) - EMAT(LLL,1)
          IF (MASWRK) THEN
            WRITE(IW,9160) III
            WRITE(IW,9180) ECONT(2)
            WRITE(IW,9190)
          END IF
CZ      ELSEIF (III.EQ.3.AND.M1M2EX(ISTEX(III-2)).GT.0) THEN
CZ        CALL POLFIT(LCUNT,2,0,XX,YY,COEF,SIGMAY,CHISQR)
CZ        ECONT(2) = COEF(1) + COEF(2)*ECONT(1)
CZ        WRITE(IW,'(/3X,A,F12.8/)') 'DELTA E(3)=', ECONT(2)
CZ      ELSEIF (III.GE.4) THEN
        ELSEIF (III.GE.3) THEN
          DELTA1 = ZERO
          DELTA2 = ZERO
          IF (ISCHME.EQ.2.AND.III.GE.5) THEN
            CALL TLMSFIT(LCUNT,XX,YY,ZZ,COEF)
            ECONT(III-1) = COEF(1)*ECONT(1) + COEF(2)*ECONT(2) + COEF(3)
            DO IIK = 1, LCUNT
              TMP1 = ZZ(IIK)-COEF(1)*XX(IIK)-COEF(2)*YY(IIK)-COEF(3)
              DELTA1 = DELTA1 + TMP1**2
              IF (ABS(TMP1).GT.DELTA2) DELTA2 = ABS(TMP1)
            END DO
          ELSE
            CALL POLFIT(LCUNT,2,0,XX,YY,COEF,SIGMAY,CHISQR)
            IF (III.EQ.3) THEN
              ECONT(III-1) = COEF(1) + COEF(2)*ECONT(1)
            ELSE
              ECONT(III-1) = COEF(1) + COEF(2)*ECONT(III-3)
            ENDIF
            DO IIK = 1, LCUNT
              TMP1 = YY(IIK)-COEF(1)-COEF(2)*XX(IIK)
              DELTA1 = DELTA1 + TMP1**2
              IF (ABS(TMP1).GT.DELTA2) DELTA2 = ABS(TMP1)
            END DO
          ENDIF
          DELTA1 = SQRT(DELTA1/LCUNT)
C
          X0 = XX(1)
          X1 = XX(LCUNT)
          X2 = ECONT(III-3)
          TMP1  = 1.0D+00 + 2.0D+00*(X2-X1)/(X1-X0)
          EMAX1 = DELTA1*TMP1
          EMAX2 = DELTA2*TMP1
C
          II1 = III
          II2 = III - 1
          IF (ISCHME.EQ.2.AND.III.GE.5) THEN
            II1 = III
            II2 = III - 2
          END IF
          IF (M1M2EX(ISTEX(III-2)).EQ.0) THEN
C
            MINT1 = MINTMP
            IF (ICUNT.GT.1) MINT1 = MINTMP - 3
            IF (DELTA1.LE.XTRTOL(III-3)) THEN
              WRITE(IW,9200) III,MINT1,MAXTMP,DELTA1,EMAX1,
     *                 DELTA2,EMAX2,LCUNT,II1,II2,ECONT(III-1)
            ELSE
              WRITE(IW,9205) III,MINT1,MAXTMP,DELTA1,EMAX1,
     *                 DELTA2,EMAX2,LCUNT,II1,II2,ECONT(III-1)
C
              MINTMP = MAXTMP + 1
              IF (IDELTM.EQ.1) THEN
                MAXTMP = MINTMP + 4
              ELSE
                MAXTMP = MINTMP + 9
              ENDIF
              IF (MAXTMP.GT.MMC) MAXTMP = MMC
              IF (MINTMP.GT.MAXTMP) GOTO 4000
              GOTO 100
            ENDIF
          ELSE
            WRITE(IW,9300) III
            WRITE(IW,9310) (M1M2EX(ISTEX(III-2)+JJ),JJ=0,MEX(III-2)-1)
            WRITE(IW,9315) LCUNT
            WRITE(IW,9320)
            IF (ISCHME.EQ.2.AND.III.GT.5) THEN
              WRITE(IW,9325) COEF(1), COEF(2), COEF(3)
            ELSE
              WRITE(IW,9330) COEF(2), COEF(1)
            END IF
            WRITE(IW,9335) DELTA1,DELTA2
            WRITE(IW,9340) II1,II2,ECONT(III-1),EMAX1,EMAX2
          ENDIF
        ENDIF
C
 4000   CONTINUE
C
CZ
CZ      WRITE OUT EMAT AND DEMAT
CZ
        IF (MASWRK) THEN
          WRITE(IW,'(/)')
          WRITE(IW,9225)
          WRITE(IW,'(5X,A)') 'CEEIS DIFFERENCE MATRIX DEMAT(M,X)'
C
          IF (ISCHME.EQ.1) THEN
            WRITE(IW,'(1X,A1,12A14)') 'M', (CDELE(I),I=1,3),
     *                             (CDELE1(I),I=1,ISTPEX-4)
          ELSE
            WRITE(IW,'(1X,A1,12A14)') 'M', (CDELE(I),I=1,3),
     *                             (CDELE2(I),I=1,ISTPEX-4)
          END IF
CZ
          NAVO = MMC - MMIN + 1
          IC = MMIN
          DO I = 1, NAVO - 1
            WRITE(IW,9215) IC, (DEMAT(I,J),J=1,ISTPEX-1)
            IC = IC + 1
          END DO
          WRITE(IW,9215) IC, (ECONT(J),J=1,ISTPEX-1)
          WRITE(IW,9225)
CZ
        ENDIF
        CALL FLSHBF(IP)
        CALL FLSHBF(IW)
      END DO
C
      ETOT = ENREF
      CTOT = ZERO
      IF (MASWRK) WRITE(IW,9350)
      DO III = 1, ISTPEX-1
        CTOT = CTOT + ECONT(III)
        IF (MASWRK) WRITE(IW,9400) III+1,ECONT(III)
      END DO
      ETOT = ETOT + CTOT
      IF (MASWRK) WRITE(IW,9450) CTOT, ENREF, ETOT
      RETURN
C
 9000 FORMAT(/3X,62(1H-)/
     *   3X, 'CALCULATIONS FOR EXCITATION LEVEL', I3,
     *   1X, 'WITH',I4,1X,'VIRTUAL ORBITALS'/
     *   3X, 'NCORE=', I4, 1X, 'NACT=', I5/
     *   3X, 62(1H-))
 9004 FORMAT(/1X,'STATE',I4,'  ENERGY= ',F20.10,'  S=',F6.2,
     *           '  SZ=',F6.2,:,'  SPACE SYM=',A4/)
 9005 FORMAT(/1X,'CEEIS ENERGY MATRIX ELEMENT EMAT(',I3,',',I2,')=',
     *       F20.10)
 9010 FORMAT(/3X,44(1H-)/
     *   3X,'  DETERMINATION OF DAVIDSON + Q CORRECTION'/
     *   3X,44(1H-))
 9020 FORMAT(/1X,'E(MR-CISD) =    ',F20.10)
 9025 FORMAT(/1X,'REF WEIGHT =    ',F20.10)
C9030 FORMAT(/1X,'*** CALCULATION OF REFERENCE ENERGY     ***')
C9040 FORMAT(/1X,'THE NUMBER OF SPACES             =',I4/
C    *        1X,'EACH SPACE STARTS AT ORBITAL     =',50I4)
 9050 FORMAT(1X,'NO OF ORBITALS IN EACH SPACE     =',50I4)
 9060 FORMAT(1X,'MIN NO OF ELECS IN EACH SPACE    =',50I4)
 9070 FORMAT(1X,'MAX NO OF ELECS IN EACH SPACE    =',50I4)
 9080 FORMAT(/1X,'E(REF)      =   ',F20.10)
 9090 FORMAT(/1X,'*** END OF REFERENCE ENERGY CALCULATION ***')
 9100 FORMAT(/1X,'E(MRCISD + Q) = ',F20.10)
 9110 FORMAT(/3X,'  END OF DAVIDSON + Q CORRECTION CALC ')
 9120 FORMAT(3X,44(1H-))
 9150 FORMAT(/1X,'THE NUMBER OF DETERMINANTS HAVING SPACE SYMMETRY ',A3/
     *        1X,'IN POINT GROUP ',A4,' WITH SZ=',F5.1,' IS',I15)
 9160 FORMAT(/3X,55(1H*)/
     *   3X,'TOTAL CONTRIBUTION FROM EXCITATION LEVEL',I2,':')
 9170 FORMAT(3X, 'DELTA E(2) = E(2) - ENREF = ', F12.8)
 9180 FORMAT(3X, 'DELTA E(3) = E(3) - E(2)  = ', F12.8)
 9190 FORMAT(3X,55(1H*)/)
 9200 FORMAT(/3X,60(1HX)/
     *   4X, 'CONVERGENCE IN THIS RANGE HAS BEEN ACHIEVED'/
     *   4X, 'EXTRAPOLATION FOR EXCITATION LEVEL',I3/
     *   4X, 'ENERGIES USED FOR LINEAR FIT M=',I4,'-',I4/
     *   4X, 'STANDARD DEVIATION OF FIT DELTA1   =',F12.8/
     *   4X, 'EXTRAPOLATION UNCERTAINTY EPSILON1 =',F12.8/
     *  /4X, 'MAXIMUM DEVIATION OF FIT DELTA2    =',F12.8/
     *   4X, 'EXTRAPOLATION UNCERTAINTY EPSILON2 =',F12.8/
     *   4X, 'NUMBER OF POINTS',I4/
     *   4X, 'EXTRAPOLATED ENERGY CONTRIBUTION'/
     *   4X,  'E(',I2,')-E(',I2,') = ',F12.8/
     *   4X,  60(1HX))
 9205 FORMAT(/3X,60(1HX)/
     *   4X, 'CONVERGENCE IN THIS RANGE HAS NOT BEEN ACHIEVED'/
     *   4X, 'EXTRAPOLATION FOR EXCITATION LEVEL',I3/
     *   4X, 'ENERGIES USED FOR LINEAR FIT M=',I4,'-',I4/
     *   4X, 'STANDARD DEVIATION OF FIT DELTA1   =',F12.8/
     *   4X, 'EXTRAPOLATION UNCERTAINTY EPSILON1 =',F12.8/
     *  /4X, 'MAXIMUM DEVIATION OF FIT DELTA2    =',F12.8/
     *   4X, 'EXTRAPOLATION UNCERTAINTY EPSILON2 =',F12.8/
     *   4X, 'NUMBER OF POINTS',I4/
     *   4X, 'EXTRAPOLATED ENERGY CONTRIBUTION'/
     *   4X,  'E(',I2,')-E(',I2,') = ',F12.8/
     *   4X,  60(1HX))
 9210 FORMAT(I2,5F16.8)
 9215 FORMAT(I2,12F14.8)
 9220 FORMAT(/1X,'FOR MXXPAN =',I4,2X,I8,' IS REQUIRED ',
     *       'BUT ACTUAL AMOUNT IS ',I8)
 9225 FORMAT(70(1H~))
 9230 FORMAT(1X,A)
 9300 FORMAT(/3X,60(1HX)/
     *   4X, 'EXTRAPOLATION FOR EXCITATION LEVEL',I3)
 9310 FORMAT(4X, 'M VALUES USED FOR LINEAR FIT M=',8I4)
 9315 FORMAT(4X, 'NUMBER OF DATA POINTS',I4)
 9320 FORMAT(/4X, 'LEAST MEAN SQUARE FIT')
 9325 FORMAT(4X, 'DE(X) =',F10.6,'DE(2) +',F10.6,'DE(3) +',F10.6)
 9330 FORMAT(4X, 'Y =',F10.6,'X +',F10.6)
 9335 FORMAT(4X, 'STANDARD DEVIATION OF FIT DELTA1 =',F12.8/
     *       4X, ' MAXIMUM DEVIATION OF FIT DELTA2 =',F12.8)
 9340 FORMAT(4X, 'EXTRAPOLATED ENERGY DIFFERENCE'/
     *   4X, 'E(',I2,')-E(',I2,') = ',F16.10/
     *  /4X, 'EXTRAPOLATION UNCERTAINTY EPSILON1 =',F12.8/
     *   4X, 'EXTRAPOLATION UNCERTAINTY EPSILON2 =',F12.8/
     *   3X, 60(1HX))
 9350 FORMAT(/3X,60(1HT)/
     *  /4X, 'TOTAL ENERGY ACCOUNTING')
 9400 FORMAT(4X, 'CONTRIBUTION FROM EXCITATION LEVEL',I3,':',F16.10)
 9450 FORMAT(4X,60(1H-)/
     *   4X, 'CONTRIBUTION OF ALL EXCITATIONS:',F16.10/
     *   4X, 'REFERENCE ENERGY               :',F16.10/
     *   4X, 'TOTAL ENERGY                   :',F16.10/
     *  /3X, 60(1HT)/)
C
      END
C*MODULE CEEIS    *DECK TLMSFIT
      SUBROUTINE TLMSFIT(NP,XX,YY,ZZ,COEF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION XX(NP), YY(NP), ZZ(NP)
      DIMENSION A(3,3), COEF(3), IPVT(3)
C
      PARAMETER (ZERO=0.0D+00)
C
CZ    LEAST MEAN SQUARE FIT (THREE VALUABLES)
CZ    ZZ = A*XX + B*YY + C
CZ    X, Y, Z: SETS OF DATA POINTS
CZ    A, B, C: VALUABLES
C
      X  = ZERO
      Y  = ZERO
      Z  = ZERO
      X2 = ZERO
      Y2 = ZERO
      XY = ZERO
      ZX = ZERO
      ZY = ZERO
      DO II = 1, NP
        X  = X  + XX(II)
        Y  = Y  + YY(II)
        Z  = Z  + ZZ(II)
        X2 = X2 + XX(II)*XX(II)
        Y2 = Y2 + YY(II)*YY(II)
        XY = XY + XX(II)*YY(II)
        ZX = ZX + ZZ(II)*XX(II)
        ZY = ZY + ZZ(II)*YY(II)
      ENDDO
C
      A(1,1)  = X2
      A(2,1)  = XY
      A(1,2)  = XY
      A(3,1)  =  X
      A(1,3)  =  X
      A(2,2)  = Y2
      A(2,3)  =  Y
      A(3,2)  =  Y
      A(3,3)  = 1.0D+00 * NP
      COEF(1) = ZX
      COEF(2) = ZY
      COEF(3) =  Z
C
      CALL SLVLEQ(A,COEF,IPVT,3,3,0,IERR)
C
      RETURN
      END
