C 20 AUG 07 - MWS - CISAO: DON'T READ DISK IF NUMBER OF ORBITALS CHANGES
C 24 MAR 06 - MWS - SKIP SINGLET/TRIPLET TRANSITIONS...THEY VANISH
C 22 DEC 06 - DGF - SYNCHRONISE MP2PAR
C  6 NOV 06 - MWS - SWITCH STATE SEPECIFIER TO USUAL KEYWORD -IROOT-
C 14 NOV 05 - SPW - CISNOS: OUTPUT CIS NATURAL ORBITALS TO LOG FILE,
C                   CISAO: FIX MEMORY PROBLEM FOR MULTI-BATCHED JOBS,
C                   CISOUT,TRDENM: USE LATEST NIST CONVERSION FACTORS
C 19 SEP 05 - MWS - FIX CIS VECTOR READING FOR SPHERICAL HARMONIC CASE
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 30 APR 05 - SPW,PA - TRDENM: ADD TRANSITION DIPOLE MOMENT CALCULATION
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL
C  7 SEP 04 - MWS - DRFCIS: MODIFY INTEGRAL STORAGE PATTERN
C 23 JUL 04 - MWS - CISINP: TIGHTEN CI CONVERGENCE IF DOING DERIVATIVES
C  7 APR 04 - MWS - CISAO: IMPOSE ORBITAL PHASES TO REUSE CI VECTOR
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  3 SEP 03 - SPW - NEW MODULE FOR CI-SINGLES ENERGY AND GRADIENT
C
C*MODULE CISGRD  *DECK CISINP
      SUBROUTINE CISINP
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL OK,PK,PANDK,BLOCK,ABEL,GOPARR,DSKWRK,MASWRK,
     *        CISPRP,MNMEDG,MNMEOP,UNVGSS,DGAPRX,RDCISV,CLOBBR
C
      PARAMETER (MXATM=2000, MXAO=8192, MXSH=5000)
C
      COMMON /CISPAR/ HAMTYP,DIAGZN,DAVCVG,PRTTOL,NSTATE,ISTATE,MULT,
     &                MXV,NDAVIT,ICISPR,NACORE,NBCORE,NOA,NOB,NORBOC,
     &                NBF,NGSVEC,MNMEDG,MNMEOP,ICLOBBR,UNVGSS,DGAPRX,
     &                RDCISV
      COMMON /CISCHF/ CHFSLV
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MP2PAR/ DUMMY(4),IDUMMY(12),MAXITC
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PKFIL / PK,PANDK,BLOCK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMMOL/ GROUP,COMPLEX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     &                TSYM(432),INVT(48),NT
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
      PARAMETER (NNAM=21)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA WRDCIS /8HCIS     /
C          NOTE: ISTATE (DEPRECATED) IS THE SAME AS IROOT KEYWORD.
      DATA QNAM   /8HNACORE  ,8HNBCORE  ,8HNSTATE  ,8HPRTTOL  ,
     &             8HISTATE  ,8HHAMTYP  ,8HMULT    ,8HDIAGZN  ,
     &             8HMXVEC   ,8HNDAVIT  ,8HDAVCVG  ,8HCISPRP  ,
     &             8HNGSVEC  ,8HMNMEDG  ,8HMNMEOP  ,8HCHFSLV  ,
     &             8HRDCISV  ,8HDGAPRX  ,8HMAXITC  ,8HCLOBBR  ,
     &             8HIROOT   /
      DATA KQNAM  /1,1,1,3,  1,5,1,5,  1,1,3,0,  1,0,0,5,
     &             0,0,1,0,  1/
C
      DATA BLANK/8H        /
      DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF    /
      DATA OPTIM,SADPT/8HOPTIMIZE,8HSADPOINT/
      DATA AIRC,DRC/8HIRC     ,8HDRC     /
      DATA C1/8HC1      /
C
      DATA SAPS,DETS/8HSAPS    ,8HDETS    /
      DATA DAVID,FULL/8HDAVID   ,8HFULL    /
      DATA CONJG,DIIS/8HCONJG   ,8HDIIS    /
C
C     --- KILL JOB IF ANYTHING OTHER THAN RHF REFERENCE ---
C
      IF(SCFTYP.EQ.ROHF .OR. SCFTYP.EQ.UHF) THEN
         IF (MASWRK) THEN
            WRITE(IW,*) 'CIS WITH UHF OR ROHF REFERENCE NOT AVAILABLE'
         END IF
         CALL ABRT
      END IF
C
C     --- SET UP FOR SCF LEVEL CIS RUNS, READ $CIS INPUT ---
C
      NORBOC= NQMT
      NBF   = NUM
      NOA   = NA
      NOB   = NB
      NVIRA = NORBOC - NOA
C
C     --- INITIALIZE NACORE (#ALPHA CORE ELECTRONS) AND
C         NBCORE (#BETA CORES) ---
C
      NACORE = NUMCOR()
      NBCORE = -1
C
C     --- SET DEFAULTS ---
C
      CISPRP = .FALSE.
      MNMEDG = .FALSE.
      MNMEOP = .FALSE.
      RDCISV = .FALSE.
      DGAPRX = .TRUE.
      CLOBBR = .FALSE.
C
      HAMTYP = BLANK
      DIAGZN = BLANK
      CHFSLV = BLANK
C
      NSTATE = 11
C          THE INPUT READS KEYWORDS ISTATE AND IROOT, TO PRESERVE
C          EXISTING INPUT FILES FROM WHEN THE SPELLING WAS ISTATE.
C          ONLY IROOT IS DOCUMENTED NOW.  THE VALUE IS CARRIED
C          THROUGH COMMON WITH THE ORIGINAL SPELLING OF ISTATE.
      ISTATE = -27
      IROOT = 1
C
      MULT = 1
      MXV = -1
      NDAVIT = 50
      DAVCVG = 1.0D-05
      CALL DERCHK(NDER)
      IF(NDER.GT.0) DAVCVG = 5.0D-06
      NGSVEC = -1
      MAXITC = 50
      PRTTOL=0.05D+00
C
      CALL NAMEIO(IR,JRET,WRDCIS,NNAM,QNAM,KQNAM,
     &            NACORE,NBCORE,NSTATE,PRTTOL,ISTATE,HAMTYP,MULT,
     &            DIAGZN,MXV,NDAVIT,DAVCVG,CISPRP,NGSVEC,MNMEDG,
     &            MNMEOP,CHFSLV,RDCISV,DGAPRX,MAXITC,CLOBBR,IROOT,
     &    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,*) 'ERROR READING $CIS GROUP'
         CALL ABRT
      END IF
      NERR = 0
C
      IF (NBCORE.EQ.-1) NBCORE = NACORE
C
C       IF ANCIENT KEYWORD ISTATE WAS NOT TYPED, USE NEW KEYWORD IROOT.
      IF(ISTATE.EQ.-27) THEN
         ISTATE=IROOT
      ELSE
         IF(MASWRK) WRITE(IW,9015)
      END IF
C
C     --- CHECK SPELLING ON CIS OPTIONS ---
C
      OK = .FALSE.
      IF(HAMTYP.EQ.BLANK) HAMTYP = SAPS
      IF(HAMTYP.EQ.SAPS)   OK=.TRUE.
      IF(HAMTYP.EQ.DETS)   OK=.TRUE.
      IF(.NOT.OK) THEN
         IF (MASWRK) WRITE(IW,9010) 'HAMTYP',HAMTYP
         NERR = NERR+1
      END IF
C
      OK = .FALSE.
      IF(DIAGZN.EQ.BLANK) DIAGZN = DAVID
      IF(DIAGZN.EQ.DAVID)   OK=.TRUE.
      IF(DIAGZN.EQ.FULL)    OK=.TRUE.
      IF(.NOT.OK) THEN
         IF (MASWRK) WRITE(IW,9010) 'DIAGZN',DIAGZN
         NERR = NERR+1
      END IF
C
      OK = .FALSE.
      IF(CHFSLV.EQ.BLANK) CHFSLV = CONJG
      IF(CHFSLV.EQ.CONJG)   OK=.TRUE.
      IF(CHFSLV.EQ.DIIS)    OK=.TRUE.
      IF(.NOT.OK) THEN
         IF (MASWRK) WRITE(IW,9010) 'CHFSLV',CHFSLV
         NERR = NERR+1
      END IF
C
C     --- SET DEFAULT FOR NUMBER OF INITIAL GUESS CIS VECTURES ---
C
      IF(NGSVEC.EQ.-1) THEN
         NGSVEC = NSTATE*2
         IF(NGSVEC.LT.10) NGSVEC = 10
      END IF
C
      NOCA  = NOA - NACORE
      IF(HAMTYP.EQ.SAPS) THEN
         NCFG = NOCA*NVIRA
      ELSE
         NCFG = 2*NOCA*NVIRA
      END IF
C
C     --- SET UP DEFAULT MAXIMUM NUMBER OF EXPANSION VECTORS ---
C
      IF(MXV.EQ.-1) THEN
         IF(NCFG.GT.1.0D+05) THEN
            MXV = 4*NSTATE
         ELSE
            MXV = 8*NSTATE
         END IF
      END IF
      IF(MXV.LT.NGSVEC) MXV = NGSVEC
C
      ICLOBBR=0
      IF(CLOBBR) ICLOBBR=1
C
      IF(DIAGZN.EQ.FULL) THEN
C
C        --- IF FULL IN MEMORY DIAGONALIZATION SET FULL SPACE ---
C
         NGSVEC = NCFG
         UNVGSS = .TRUE.
         RDCISV = .FALSE.
C
      ELSE IF(RDCISV) THEN
C        INITIAL CIS VECTORS READ FROM INPUT FILE IN WFNCIS
         UNVGSS = .FALSE.
      ELSE
C        UNIT VECTOR GUESS
         UNVGSS = .TRUE.
      END IF
C
C     --- OPTIONAL PROPERTIES DURING SINGLE POINT ENERGY RUN ---
C
      IF(CISPRP) THEN
         ICISPR = 1
      ELSE
         ICISPR = 0
      END IF
C
      IF(MASWRK) WRITE(IW,9020) NACORE,NBCORE,NSTATE,ISTATE,
     &                          HAMTYP,MULT,DIAGZN,MXV,NDAVIT,
     &                          DAVCVG,CISPRP,NGSVEC,MNMEDG,
     &                          MNMEOP,CHFSLV,RDCISV,DGAPRX
C
C     --- FOR NOW IF GEOMETRY SEARCH SWITCH OFF SYMMETRY ---
C
      IF(NDER.GT.0) THEN
         ICISPR = 0
         IF(RUNTYP.EQ.OPTIM .OR. RUNTYP.EQ.SADPT .OR.
     &      RUNTYP.EQ.AIRC  .OR. RUNTYP.EQ.DRC) THEN
            GROUP = C1
            IGROUP = 1
            NAXIS = 1
            NT = 1
            IF(MASWRK) WRITE(IW,9025)
         END IF
      END IF
C
C     --- MORE INCOMPATABILITY TESTS ---
C
      IF(SCFTYP.NE.RHF) THEN
         IF(MASWRK) WRITE(IW,9090)
         NERR = NERR+1
      END IF
      IF(GZRF.NE.0.0D+00) THEN
         IF(MASWRK) WRITE(IW,9120)
         NERR = NERR+1
      END IF
C
C     --- FILL UP PAIR INDEX ARRAY IN COMMON ---
C
      DO 5 I=1,NBF
         IA(I) = (I*I-I)/2
    5 CONTINUE
C
      IF(MASWRK) WRITE(IW,9030) NACORE,NBCORE,NOA,NOB,NORBOC,NBF
C
      IF(NORBOC.LE.NOA) THEN
         IF (MASWRK) WRITE(IW,9040)
         NERR = NERR+1
      END IF
      IF(NORBOC.LE.NACORE .OR. NORBOC.LE.NBCORE) THEN
         IF(MASWRK) WRITE(IW,9060) NACORE,NBCORE,NORBOC
         NERR = NERR+1
      END IF
      IF(NACORE.GT.NOA  .OR.  NBCORE.GT.NOB .OR. (NACORE.EQ.NOA
     &   .AND.NBCORE.EQ.NOB)) THEN
         IF(MASWRK) WRITE(IW,9065) NACORE,NBCORE,NOA,NOB
         NERR = NERR+1
      END IF
C
      NOPK = 1
      PK   = .FALSE.
      PANDK= .FALSE.
C
      IF(NERR.GT.0) THEN
         IF(MASWRK) WRITE(IW,9050)
         CALL ABRT
      END IF
C
      RETURN
C
 9010 FORMAT(/1X,'ERROR: $CIS KEYWORD ',A6,
     &          ' WAS GIVEN AN ILLEGAL VALUE ',A8,'.')
 9015 FORMAT(/1X,'*** WARNING ***'/
     *        1X,'OBSOLETE KEYWORD SPELLING -ISTATE- FOUND IN $CIS,'/
     *        1X,'FROM NOW ON, PLEASE GIVE SPELL THIS AS -IROOT-'/)
 9020 FORMAT(/5X,30(1H-)/5X,'CI-SINGLES CONTROL INFORMATION'/5X,30(1H-)/
     &        5X,'NACORE = ',I8,2X,'NBCORE   = ',I8/
     &        5X,'NSTATE = ',I8,2X,'IROOT    = ',I8/
     &        5X,'HAMTYP = ',A8,2X,'SAP MULT = ',I8/
     &        5X,'DIAGZN = ',A8,2X,'MXVEC    = ',I8/
     &        5X,'NDAVIT = ',I8,2X,'DAVCVG   = ',1P,E8.2/
     &        5X,'CISPRP = ',L8,2X,'NGSVEC   = ',I8/
     &        5X,'MNMEDG = ',L8,2X,'MNMEOP   = ',L8/
     &        5X,'CHFSLV = ',A8,2X,'RDCISV   = ',L8/
     &        5X,'DGAPRX = ',L8)
 9025 FORMAT(/5X,'SWITCHING OFF SYMMETRY FOR CIS GEOMETRY SEARCH')
 9030 FORMAT(/5X,'NUMBER OF CORE -A-  ORBITALS = ',I5/
     &        5X,'NUMBER OF CORE -B-  ORBITALS = ',I5/
     &        5X,'NUMBER OF OCC. -A-  ORBITALS = ',I5/
     &        5X,'NUMBER OF OCC. -B-  ORBITALS = ',I5/
     &        5X,'NUMBER OF MOLECULAR ORBITALS = ',I5/
     &        5X,'NUMBER OF   BASIS  FUNCTIONS = ',I5/)
 9040 FORMAT(1X,'***ERROR*** CIS HAS NO VIRTUAL ORBITALS.')
 9050 FORMAT(5X,'*** PROGRAM STOPPED IN -CISINP- ***')
 9060 FORMAT(' INCONSISTENT INPUT PARAMETERS IN -$CIS- .',
     &       ' NACORE,NBCORE,NORBOC = ',3I5)
 9065 FORMAT(1X,'*** ERROR DETECTED ***'/
     &       1X,'NACORE=',I5,' NBCORE=',I5,' OCCUPIED ALPHA,BETA=',2I5//
     &       1X,'IF THE FROZEN CORE CONTAINS ALL ELECTRONS, OR'/
     &       1X,'IF NACORE OR NBCORE IS GREATER THAN THE NUMBER OF'/
     &       1X,'ALPHA OR BETA ELECTRONS,'/
     &       1X,'PLEASE RETHINK YOUR RUN PARAMETERS'//)
 9090 FORMAT(1X,'*** ERROR *** ',
     &          'CIS ENERGY IS AVAILABLE ONLY FOR RHF REFERENCE')
 9120 FORMAT(1X,'*** ERROR *** ',
     &          'CIS IS PRESENTLY INCOMPATIBLE WITH USE OF SCRF')
      END
C
C*MODULE CISGRD  *DECK CISAO
      SUBROUTINE CISAO
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30,
     &           MXRT=100, MXATM=2000, MXAO=8192)
C
      LOGICAL PACK2E,DIRSCF,FDIFF,GOPARR,DSKWRK,MASWRK,
     &        ABEL,ABELPT,DETRUN,UNVGSS,DGAPRX,RDCISV,
     &        MNMEDG,MNMEOP,CISSNG,TDSKWRK
C
      DIMENSION NOCP(1)
C
      COMMON /CIFILS/ NFT11,NFT12,NFT13,NFT14,NFT15,NFT16,IDAF20,NEMEMX
      COMMON /CISCHF/ CHFSLV
      COMMON /CISPAR/ HAMTYP,DIAGZN,DAVCVG,PRTTOL,NSTATE,ISTATE,MULT,
     &                MXV,NDAVIT,ICISPR,NACORE,NBCORE,NOA,NOB,NORBOC,
     &                NBF,NGSVEC,MNMEDG,MNMEOP,ICLOBBR,UNVGSS,DGAPRX,
     &                RDCISV
      COMMON /CISEN / ECIS
      COMMON /CISMLT/ CISSNG
      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
      COMMON /DIRMEM/ IPQKL,IVEC,IDUM1,IDUM2,NUM,NOCC,NDIM,MINPQ,MAXPQ
      COMMON /ENRGYS/ ENUC,EELEC,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     &                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FMCOM / X(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NNUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     &                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     &                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     &                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMMOL/ GROUP,COMPLEX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     &                TSYM(432),INVT(48),NT
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00,
     *           DEBYE=2.541766D+00)
C
      DATA SAPS,DETS/8HSAPS    ,8HDETS    /
      DATA CONJG,DIIS/8HCONJG   ,8HDIIS    /
      DATA CHECK/8HCHECK   /
C
      SAVE NCFGSAVE
C
C     DRIVER FOR ATOMIC-ORBITAL BASIS CI-SINGLES ENERGY AND GRADIENT
C         WRITTEN BY SIMON P. WEBB, JUNE 2003.
C
      CALL DERCHK(NDER)
C
      NOCP(1) = 0
C
C     --- LATER CALCULATION OF TRANSITION DIPOLE MOMENT REQUIRES
C     DIPOLE INTEGRALS AT CURRENT GEOMETRY, WHICH MAY OR MAY NOT HAVE
C     BEEN COMPUTED EARLIER IN THE RUN, AND NOT NECESSARILY AT A KNOWN
C     ORIGIN. DO THIS NOW. ---
C
      CALL CALCOM(XP,YP,ZP)
      CALL DIPINT(XP,YP,ZP,.FALSE.)
C
C     --- ORBITAL AND BASIS FUNCTION PARAMETERS ---
C
      NOCA  = NOA - NACORE
CSPW  NVIRA = NORBOC - NOA
      NOC   = NOA
CSPW  NVIR  = NVIRA
      NVIR  = NQMT - NOC
C
      NBF2  = (NBF*NBF+NBF)/2
      NBF3  = NBF*NBF
      NUM   = NBF
      NDIM  = NBF2
C
C     --- GET CORRECT INTEGRAL SCALE FACTOR FOR DETS OR SAPS ---
C
      CISSNG = .TRUE.
      DETRUN = .FALSE.
      NALBE = 1
C
      IF(HAMTYP.EQ.DETS) THEN
         SCLINT = ONE
         DETRUN = .TRUE.
         NALBE = 2
      ELSE
         IF(MULT.EQ.1) THEN
            SCLINT = TWO
         ELSE IF(MULT.EQ.3) THEN
            SCLINT = ZERO
            CISSNG = .FALSE.
         ELSE
            IF(MASWRK) WRITE(IW,*) 'PROBLEM WITH MULT IN CISAO - ABORT'
            CALL ABRT
         END IF
      END IF
C
C     --- NUMBER OF CONFIGURATIONS DEPENDS ON HAMILTONIAN TYPE ---
C
      IF(HAMTYP.EQ.SAPS) THEN
         NCFG = NOCA*NVIR
      ELSE
         NCFG = 2*NOCA*NVIR
      END IF
      NROTA = NOCA*NVIR
      NROT = NOC*NVIR
C
      IF(NSTATE.EQ.NCFG) THEN
         IF(MASWRK) WRITE(IW,9010)
         MXV = NCFG
         NGSVEC = NCFG
      END IF
      IF(MXV.GT.NCFG) THEN
         IF(MASWRK) WRITE(IW,9015)
         MXV = NCFG
         NGSVEC = NCFG
      END IF
C
      MXV2 = (MXV*MXV+MXV)/2
C
      IF(.NOT.UNVGSS) NGSVEC = NSTATE
C
      IF(MASWRK) THEN
         IF(NDER.EQ.0  .AND.  ICISPR.EQ.0) THEN
            WRITE(IW,9000)
         ELSE IF(NDER.EQ.0  .AND.  ICISPR.EQ.1) THEN
            WRITE(IW,9001)
         ELSE
            WRITE(IW,9002)
         END IF
         WRITE(IW,9005) NACORE,NOCA,NQMT,NBF
         IF(HAMTYP.EQ.SAPS) THEN
            WRITE(IW,9006) NCFG
         ELSE
            WRITE(IW,9007) NCFG
         END IF
      END IF
C
C     --- GET MEMORY TO BE USED THROUGHOUT ---
C
      CALL VALFM(LOADFM)
      LVEC  = 1     + LOADFM
      IIRP  = LVEC  + NBF3
      IENG  = IIRP  + NBF
      LAST  = IENG  + NBF
      NEEDA = LAST  - LOADFM
      CALL GETFM(NEEDA)
C
      IF(EXETYP.EQ.CHECK) GO TO 90
C
C     --- READ IN ORBITALS AND ORBITAL ENERGIES ---
C     A STANDARD PHASE IS IMPOSED ON THE ORBITALS, TO FACILITATE
C     THE RECYCLING OF CIS VECTORS FROM PREVIOUS GEOMETRIES.
C
      CALL DAREAD(IDAF,IODA,X(LVEC),NBF3,15,0)
      CALL STFASE(X(LVEC),NBF,NBF,NQMT)
      CALL DAREAD(IDAF,IODA,X(IENG),NBF,17,0)
C
   90 CONTINUE
C
CSPW  IRREPS IN X(IIRP) MAY USED IN THE FUTURE
CSPW
CSPW  --- READ IN OVERLAP INTEGRALS, Q-MATRIX, AND DETERMINE IREPS ---
CSPW
CSPW  CALL VALFM(LOADFM)
CSPW  IWRK   = 1     + LOADFM
CSPW  IWRK2  = IWRK  + NBF
CSPW  IOVLP  = IWRK2 + NBF
CSPW  IQMAT  = IOVLP + NBF2
CSPW  ITMP   = IQMAT + NBF3
CSPW  LAST   = ITMP  + NBF
CSPW  NEEDB  = LAST  - LOADFM
CSPW  CALL GETFM(NEEDB)
CSPW
CSPW  CALL DAREAD(IDAF,IODA,X(IOVLP),NBF2,12,0)
CSPW  CALL DAREAD(IDAF,IODA,X(IQMAT),NBF3,45,0)
CSPW
CSPW  CALL TRFSYM(X(IWRK),X(IIRP),X(IWRK2),X(IQMAT),X(IOVLP),
CSPW &            X(LVEC),X(ITMP),IA,NQMT,NBF,NORBOC,NBF)
CSPW
CSPW  CALL RETFM(NEEDB)
C
C     --- MEMORY ALLOCATION FOR DIRECT CI-SINGLES CALCULATION ---
C
      CALL VALFM(LOADFM)
      LWRK1  = 1      + LOADFM
      LWRK2  = LWRK1  + NBF3
      NBUF   = LWRK2  + NBF3
      IBUF   = NBUF   + NINTMX
      LAI    = IBUF   + NINTMX
      LAST   = LAI    + MAX(NCFG,NVIR*NOC)
      NEEDC  = LAST   - LOADFM
      CALL GETFM(NEEDC)
C
      ABEL = ABELPT()
      NTSAVE = NT
C
      IF(DGAPRX) THEN
C
C        --- IF APPROXIMATE DIAGONAL CIS HAMILTONIAN ELEMENTS USED
C            NO EXTRA MEMORY NEEDED ---
C
         IF(MASWRK) THEN
            WRITE(IW,9295)
         END IF
         NVBTCH = 0
         NBATCH = 1
C
      ELSE
C
C        --- IF NONABELIAN SWITCH OFF SYMMETRY ---
C
         IF(.NOT.ABEL) THEN
            IF(MASWRK) WRITE(IW,9297)
            NT = 1
            IF(.NOT.DIRSCF) CALL JANDK
         END IF
C
C        --- CHECK HOW MUCH MEMORY AVAILABLE AND SET UP BATCHING ---
C
         IF(MNMEDG) THEN
            NGOTMX = 2*NBF2
         ELSE
            CALL GOTFM(NGOTMX)
         END IF
         MEMREQ = 2*NBF2*NOCA
         IF(MEMREQ.LT.NGOTMX) THEN
            NBATCH = 1
            NVBTCH = NOCA
         ELSE
            NBATCH = MEMREQ/NGOTMX
            NVBTCH = NGOTMX/(2*NBF2)
            IF(MOD(MEMREQ,NGOTMX).NE.0) NBATCH = NBATCH + 1
         END IF
C
         IF(MASWRK) THEN
            WRITE(IW,9300)
            WRITE(IW,9310) 2*NBF2
            WRITE(IW,9320) 2*NBF2*NOCA
            WRITE(IW,9330) NGOTMX
            IF(MEMREQ.LT.NGOTMX) THEN
               WRITE(IW,9340)
            ELSE
               WRITE(IW,9350) NBATCH
            END IF
         END IF
C
      END IF
C
      CALL VALFM(LOADFM)
      LPMN   = 1      + LOADFM
      LFLM   = LPMN   + NBF2*NVBTCH
      LAST   = LFLM   + NBF2*NVBTCH
      NEEDC2 = LAST   - LOADFM
      CALL GETFM(NEEDC2)
C
      CALL TSECND(TIM0)
C
      IF(EXETYP.EQ.CHECK) GO TO 180
C
C     --- CALCULATE DIAGONAL ELEMENTS OF CI-SINGLES HAMILTONIAN ---
C
      CALL CISDGL(X(LAI),X(LVEC),X(LWRK1),X(LWRK2),X(LPMN),X(LFLM),
     &            X(NBUF),X(IBUF),X(IENG),NINTMX,
     &            NBF,NBF2,NOC,NVIR,SCLINT,NCFG,NACORE,NOCA,IA,
     &            NBATCH,NVBTCH,DGAPRX)
      IF(DETRUN) CALL DCOPY(NROTA,X(LAI),1,X(LAI+NROTA),1)
C
  180 CONTINUE
C
      NT = NTSAVE
C
      CALL RETFM(NEEDC2)
C
      CALL TSECND(TIM1)
      TIME = TIM1 - TIM0
      IF(MASWRK) WRITE(IW,9200) TIME
C
C     --- MEMORY ALLOCATION FOR DIRECT CI-SINGLES CALCULATION ---
C
      CALL VALFM(LOADFM)
      LWRK3  = 1      + LOADFM
      LPMN   = LWRK3  + NBF3
      LPIJ   = LPMN   + NBF2
      LPAB   = LPIJ   + NOC*NOC
      LFOCC  = LPAB   + NVIR*NVIR
      LFVIR  = LFOCC  + NOC*NOC
      LC1    = LFVIR  + NVIR*NVIR
      LC2    = LC1    + NVIR*NOC
      LS1    = LC2    + NVIR*NOC
      LS2    = LS1    + NOC*NOC
      LDINTX = LS2    + NVIR*NVIR
      LDINTY = LDINTX + NBF2
      LDINTZ = LDINTY + NBF2
      LAST   = LDINTZ + NBF2
      NEEDD  = LAST   - LOADFM
      CALL GETFM(NEEDD)
C
C     --- REUSE SOME MEMORY ---
C
      LWIJ   = LFOCC
      LWAB   = LFVIR
      LWAI   = LC1
C
      CALL VALFM(LOADFM)
      LA     = 1      + LOADFM
      LB     = LA     + NCFG*MXV
      LH1    = LB     + NCFG*MXV
      LH2    = LH1    + MXV2
      LDVEV  = LH2    + MXV2
      LDVVEC = LDVEV  + MXV
      LSCR1  = LDVVEC + MXV*MXV
      LSCR2  = LSCR1  + MXV*8
      LSCR3  = LSCR2  + MXV
      LRNORM = LSCR3  + MXV
      LSPN   = LRNORM + NSTATE
      LAST   = LSPN   + NSTATE
      NEEDE  = LAST   - LOADFM
      CALL GETFM(NEEDE)
C
C     --- CHECK HOW MUCH MEMORY AVAILABLE AND SET UP BATCHING ---
C         NOTE: NBATCH IS UPPER LIMIT BASED ON FACT NGSVEC.GE.NSTATE
C
      IF(MNMEOP) THEN
         NGOTMX = NBF3*NALBE*2
      ELSE
         CALL GOTFM(NGOTMX)
      END IF
      MEMREQ = NGSVEC*NBF3*NALBE*2
      IF(MEMREQ.LT.NGOTMX) THEN
         NBATCH = 1
         NVBTCH = NGSVEC
      ELSE
         NBATCH = MEMREQ/NGOTMX
         NVBTCH = NGOTMX/(NBF3*NALBE*2)
         IF(MOD(MEMREQ,NGOTMX).NE.0) NBATCH = NBATCH + 1
      END IF
C
      IF(MASWRK) THEN
         WRITE(IW,9400) NBF3*NALBE*2
         WRITE(IW,9410) NGSVEC*NBF3*NALBE*2
         WRITE(IW,9420) NGOTMX
         IF(MEMREQ.LT.NGOTMX) THEN
            WRITE(IW,9430)
         ELSE
            WRITE(IW,9440) NBATCH
         END IF
      END IF
C
      CALL VALFM(LOADFM)
      LTA    = 1      + LOADFM
      LFLA   = LTA    + NBF3*NVBTCH
      LAST   = LFLA   + NBF3*NVBTCH
      IF(DETRUN) THEN
         LTB  = LAST
         LFLB = LTB   + NBF3*NVBTCH
         LAST = LFLB  + NBF3*NVBTCH
      ELSE
         LTB  = LAST
         LFLB = LAST
      END IF
      NEEDE2 = LAST   - LOADFM
      CALL GETFM(NEEDE2)
C
C     --- NOTE TO PROGRAMMERS: UNTIL NEEDE2 MEMORY IS RETURNED DO NOT
C         ALLOCATE ANY FURTHER MEMORY - A MULTI BATCH JOB WILL ALREADY
C         HAVE TAKEN IT ALL!
C
      CALL TSECND(TIM0)
      IF(EXETYP.EQ.CHECK) GO TO 270
C
C     --- POSSIBLY READ GUESS CIS VECTORS FROM DISK FILE 12 ---
C     AFTER INITIAL GEOMETRY, THE BEST GUESS IS READING THE PREVIOUS
C     GEOMETRY'S CI VECTOR.  THE NUMBER OF CI COEFFICIENTS MIGHT CHANGE
C     IF THERE IS A CHANGE IN HOW MANY ORBITALS ARE DROPPED FOR LINEAR
C     DEPENDENCY.  THESE MUST NOT TRY TO READ PREVIOUS CI VECTOR!
C
      IF(NEVALS.EQ.0) NCFGSAVE=0
      IF(NEVALS.GT.0  .AND.  NCFG.NE.NCFGSAVE) UNVGSS=.TRUE.
      NCFGSAVE = NCFG
C
      IF(.NOT.UNVGSS) THEN
         TDSKWRK = DSKWRK
         DSKWRK = .FALSE.
         CALL SEQOPN(NFT12,'CIVECTR','UNKNOWN',.FALSE.,'UNFORMATTED')
         CALL SEQREW(NFT12)
C
         DO IST=1,NSTATE
            LPOINT = LA + (IST-1)*NCFG
            CALL SQREAD(NFT12,X(LPOINT),NCFG)
            IF(NCFG.EQ.0) THEN
               IF(MASWRK) WRITE(IW,*)
     &            'UNEXPECTED END OF FILE ON UNIT',NFT12
               CALL ABRT
            END IF
         END DO
         DSKWRK = TDSKWRK
      END IF
C
C     --- FOR NOW ALL SYMMETRY SWITCHED OFF FROM THIS POINT ON ---
C
      IF(NT.NE.1) THEN
         NT = 1
         IF(.NOT.DIRSCF .AND. (ABEL.OR.DGAPRX)) CALL JANDK
      END IF
C
C     --- CALCULATE CI-SINGLES VECTORS AND ENERGY ---
C
      CALL CISAOE(X(LVEC),X(IENG),X(NBUF),X(IBUF),
     &            X(LAI),X(LWRK1),X(LWRK2),
     &            X(LA),X(LB),X(LH1),X(LH2),X(LDVEV),X(LDVVEC),
     &            X(LSCR1),X(LSCR2),X(LSCR3),X(LRNORM),
     &            X(LTA),X(LTB),X(LFLA),X(LFLB),
     &            ESCF,NBF,NBF3,NOC,NVIR,MXV,MXV2,NINTMX,
     &            NSTATE,ISTATE,DAVCVG,NDAVIT,IA,ETOT,SCLINT,
     &            NCFG,DETRUN,NACORE,NOCA,NGSVEC,UNVGSS,NBATCH,
     &            NVBTCH,DGAPRX,IW,IP,MASWRK)
C
      ECIS = ETOT
      DO IST=1,MIN(NSTATE,MXRT)
         ESTATE(IST) = ESCF + X(LDVEV-1+IST)
      ENDDO
C
C     --- WRITE CIS VECTORS TO DISK FILE 12 ---
C
      TDSKWRK = DSKWRK
      DSKWRK = .FALSE.
      CALL SEQOPN(NFT12,'CIVECTR','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQREW(NFT12)
C
      DO  IST=1,NSTATE
         LPOINT = LA + (IST-1)*NCFG
CSPW     CALL STFASE(X(LPOINT),NCFG,NCFG,1)
         CALL SQWRIT(NFT12,X(LPOINT),NCFG)
      END DO
      CALL SEQREW(NFT12)
      DSKWRK = TDSKWRK
C
C        OPTION TO LET NEXT GEOMETRY READ CI VECTORS FROM DISK,
C        RATHER THAN UNIT VECTOR GUESS PROCEDURE.
C
      UNVGSS = .FALSE.
      IF(ICLOBBR.EQ.1) UNVGSS = .TRUE.
C
C     --- PUNCH CIS VECTORS TO .DAT FILE ---
C
      IF(MASWRK) CALL CISVPU(IP,DETRUN,NSTATE,NCFG,ESTATE,X(LA))
C
C     --- IF DETERMINANT RUN WORK OUT SPIN OF EACH STATE;
C     IF SAPS RUN JUST SET ACCORDING TO MULT ---
C
      CALL CISSPN(X(LA),X(LSPN),NROTA,NCFG,NSTATE,DETRUN,MULT)
C
C     --- OUTPUT ENERGIES AND LARGEST CIS COEFFICENTS ---
C
      IF(MASWRK) THEN
         CALL CISOUT(X(LA),X(LDVEV),X(LSPN),X(LWRK1),ESCF,NROTA,NCFG,
     &               NSTATE,NACORE,NOCA,NOC,NVIR,PRTTOL,DETRUN,IW)
      END IF
C
      IOFST = (ISTATE-1)*NCFG
C
      CALL TSECND(TIM2)
      TIME = TIM2 - TIM1
      IF(MASWRK) WRITE(IW,9205) TIME
C
C     --- FORM OCC-OCC AND VIRT-VIRT BLOCKS OF THE CI-SINGLES
C         DIFFERENCE DENSITY MATRICES FOR CHOSEN EXCITED STATE ---
C
      CALL DGEMM('T','N',NOCA,NOCA,NVIR,-ONE,X(LA+IOFST),NVIR,
     &           X(LA+IOFST),NVIR,ZERO,X(LPIJ),NOCA)
      CALL MRARTR(X(LA+IOFST),NVIR,NVIR,NOCA,X(LA+IOFST),NVIR,
     &            NVIR,X(LPAB),NVIR)
C
      IF(DETRUN) THEN
C
C        --- FOR DETERMINANT RUN ADD BETA DENSITY TO ALPHA ---
C
         CALL DGEMM('T','N',NOCA,NOCA,NVIR,-ONE,X(LA+IOFST+NROTA),NVIR,
     &              X(LA+IOFST+NROTA),NVIR,ONE,X(LPIJ),NOCA)
         CALL DGEMM('N','T',NVIR,NVIR,NOCA,ONE,X(LA+IOFST+NROTA),NVIR,
     &              X(LA+IOFST+NROTA),NVIR,ONE,X(LPAB),NVIR)
      END IF
C
C     --- IF DETERMINANT RUN TAKE APPROPRIATE LINEAR COMBINATIONS OF
C         CIS COEFFICIENTS FOR TARGET STATE SPIN AND RESET SCLINT ---
C
      IF(DETRUN) THEN
         CALL CISCMB(X(LA),X(LSPN),NROTA,NCFG,NSTATE,ISTATE,SCLINT,
     &               CISSNG)
      END IF
C
C     --- FORM TRANSITION DENSITY (X(LTA)), FOCK-LIKE MATRIX (X(LFLA)),
C         AND B MATRIX (X(LB)) FOR CONVERGED ISTATE ---
C
      CALL TFBAO(X(LVEC),X(NBUF),X(IBUF),X(LWRK1),X(LWRK2),
     &           X(LA+IOFST),X(LB),X(LTA),X(LTB),X(LFLA),X(LFLB),
     &           NCFG,NBF,NBF3,NOC,NVIR,NROTA,NINTMX,SCLINT,
     &           .FALSE.,NACORE,NOCA,0,NOC,1)
      CALL DAWRIT(IDAF,IODA,X(LTA),NBF3,301,0)
C
  270 CONTINUE
      CALL RETFM(NEEDE2)
      IF(EXETYP.EQ.CHECK) GO TO 275
C
C     --- GET DIPOLE INTEGRALS OVER AOS, AT THE CENTER OF MASS, AND THEN
C     EVALUATE THE DIPOLE MOMENT OF THE GROUND STATE DETERMINANT. ---
C
      CALL DAREAD(IDAF,IODA,X(LDINTX),NBF2,95,0)
      CALL DAREAD(IDAF,IODA,X(LDINTY),NBF2,96,0)
      CALL DAREAD(IDAF,IODA,X(LDINTZ),NBF2,97,0)
C         RHF LEVEL DENSITY MATRIX AT -WRK1-
      CALL DAREAD(IDAF,IODA,X(LWRK1), NBF2,16,0)
C         ELECTRONIC PART OF DIPOLE
      DMX = -TRACEP(X(LWRK1),X(LDINTX),NBF)
      DMY = -TRACEP(X(LWRK1),X(LDINTY),NBF)
      DMZ = -TRACEP(X(LWRK1),X(LDINTZ),NBF)
C         NUCLEAR PART OF DIPOLE
      DO IAT=1,NAT
         DMX = DMX + ZAN(IAT)*(C(1,IAT)-XP)
         DMY = DMY + ZAN(IAT)*(C(2,IAT)-YP)
         DMZ = DMZ + ZAN(IAT)*(C(3,IAT)-ZP)
      ENDDO
      DMX = DMX * DEBYE
      DMY = DMY * DEBYE
      DMZ = DMZ * DEBYE
C
C     --- CALCULATE TRANSITION DIPOLE MOMENTS ---
C
      CALL TRDENM(X(LVEC),X(LWRK1),X(LA),X(LWRK2),X(LS1),X(LS2),
     &            X(LDINTX),X(LDINTY),X(LDINTZ),X(LDVEV),X(LSPN),
     &            NCFG,NBF,NBF3,NOC,NVIR,NACORE,NOCA,NSTATE,
     &            ESCF,MASWRK,IW,DMX,DMY,DMZ)
C
C     --- FORM INTERMEDIATE MATRICES REQUIRED FOR CI-SINGLES
C         GRADIENT OR PROPERTIES CALCULATION (C1,C2,S1,S2) ---
C
      IF(NACORE.GT.0) THEN
         CALL DCOPY(NOCA*NOCA,X(LPIJ),1,X(LWRK3),1)
         CALL VCLR(X(LPIJ),1,NOC*NOC)
      END IF
C
      CALL CISITM(X(LVEC),X(LFLA),X(LWRK2),X(LFOCC),X(LFVIR),
     &            X(LA+IOFST),X(LB),X(LC1),X(LC2),X(LS1),X(LS2),
     &            X(LPIJ),X(LWRK3),X(IENG),NOC,NVIR,NBF,NACORE,NOCA)
C
  275 CONTINUE
C
      CALL RETFM(NEEDE)
C
C     --- CHECK FOR GRADIENT OR PROPERTY RUN - IF NONE TRANSFORM
C         EXPECTATION VALUE DENSITY AND WRITE TO DAF ---
C
      IF(NDER.EQ.0 .AND. ICISPR.EQ.0) THEN
C
         CALL CISAOD(X(LPIJ),X(LPAB),X(LAI),X(LVEC),X(LWRK1),X(LWRK2),
     &               X(LPMN),NBF,NBF2,NOC,NVIR,.FALSE.)
         CALL DAWRIT(IDAF,IODA,X(LPMN),NBF2,307,0)
C
C        --- ADD EXP VAL CI-SINGLES AO DIFFERENCE DENSITY TO HF DENSITY
C
         CALL DAREAD(IDAF,IODA,X(LWRK1),NBF2,16,0)
         CALL VADD(X(LPMN),1,X(LWRK1),1,X(LWRK2),1,NBF2)
         CALL DAWRIT(IDAF,IODA,X(LWRK1),NBF2,308,0)
         CALL DAWRIT(IDAF,IODA,X(LWRK2),NBF2,16,0)
C
C        --- SKIP RESPONSE DENSITY AND OTHER GRADIENT RELATED STUFF ---
C
         GO TO 540
C
      END IF
C
      IF(EXETYP.EQ.CHECK) GO TO 360
C
C     --- FORM CI-SINGLES LAGRANGAIN MATRIX ---
C
      CALL CISAOL(X(LAI),X(LPIJ),X(LPAB),X(LPMN),X(LVEC),
     &            X(LWRK1),X(LWRK2),X(LWRK3),X(NBUF),X(IBUF),
     &            X(LC1),X(LC2),NOC,NVIR,NBF,NINTMX)
C
  360 CONTINUE
C
      CALL TSECND(TIM3)
      TIME = TIM3 - TIM2
      IF(MASWRK) WRITE(IW,9210) TIME
C
C     --- SOLVE COUPLED-PERTURBED HARTREE-FOCK EQS. IN THE AO BASIS ---
C
      CALL VALFM(LOADFM)
      IPAI  = 1     + LOADFM
      IPBJ  = IPAI  + NROT
      IERR  = IPBJ  + NROT
      LAST  = IERR  + NROT
C
C     --- SET UP MEMORY FOR CHOSEN SOLVER ---
C
      IF(CHFSLV.EQ.CONJG) THEN
         IF(DIRSCF) THEN
            CALL BASCHK(LMAX)
                          NANGM =  4
            IF(LMAX.EQ.2) NANGM =  6
            IF(LMAX.EQ.3) NANGM = 10
            IF(LMAX.EQ.4) NANGM = 15
            IF(LMAX.EQ.5) NANGM = 21
            IF(LMAX.EQ.6) NANGM = 28
            MAXG = NANGM**4
            NSH2 = (NSHELL*NSHELL+NSHELL)/2
            MXG2 = MXGSH*MXGSH
         ELSE
            MAXG = 1
            NSH2 = 1
            MXG2 = 1
         END IF
C
         LPRE   = LAST
         LGHOND = LPRE   + NROT
         LXINTS = LGHOND + MAXG
         LDSH   = LXINTS + NSH2
         LDDIJ  = LDSH   + NSH2
         LAST   = LDDIJ  + 49*MXG2
      ELSE
         LPRE   = LAST
         LGHOND = LAST
         LXINTS = LAST
         LDSH   = LAST
         LDDIJ  = LAST
      END IF
C
      IF(CHFSLV.EQ.DIIS) THEN
         MAXIO   = 2*MAXDII
         MAXIT   = 50
         ITERMX  = 4*MAXIT
         ITERMX2 = (ITERMX*ITERMX+ITERMX)/2
C
         IAM    = LAST
         IX2    = IAM    + MAXDII*MAXDII
         IPVT   = IX2    + ITERMX
         IB     = IPVT   + ITERMX
         IODIIS = IB     + ITERMX2
         LAST   = IODIIS + MAXIO
      ELSE
         IAM    = LAST
         IX2    = LAST
         IPVT   = LAST
         IB     = LAST
         IODIIS = LAST
      END IF
C
      NEEDF = LAST  - LOADFM
      CALL GETFM(NEEDF)
C
      IF(MASWRK) WRITE(IW,9100)
C
      IF(EXETYP.EQ.CHECK) GO TO 450
C
      IF(CHFSLV.EQ.CONJG) THEN
C
C        --- CONJUGATE GRADIENT CPHF SOLVER (DEFAULT FOR CIS) ---
C
         NNXYZ = 1
         NUNIQ = 1
         NFOCK = 0
C
         CALL TRPOSE(X(LAI),X(LWRK1),NVIR,NOC,1)
         CALL AOCPCG(X(LAI),X(IPAI),X(IPBJ),X(IERR),X(LC1),
     &               X(LPRE),BNORM,BKNUM,BKDEN,
     &               X(LVEC),X(IENG),X(LPMN),X(LWRK1),
     &               X(LWRK2),X(LWRK3),NOCP,X(LGHOND),X(LXINTS),
     &               X(LDSH),X(LDDIJ),X(NBUF),X(IBUF),
     &               NNXYZ,NUNIQ,NFOCK,NROT,NOC,NVIR,
     &               NBF,NBF2,NBF3,NSH2,MAXG,MXG2)
C
         CALL DSCAL(NROT,-FOUR,X(IPAI),1)
         CALL TRPOSE(X(IPAI),X(LWRK1),NOC,NVIR,1)
         CALL DCOPY(NROT,X(IPAI),1,X(LAI),1)
C
      ELSE
C
C        --- DIIS CPHF SOLVER ---
C
         CALL AOCPHF(X(LPMN),X(LWRK1),X(LWRK2),X(LWRK3),X(NBUF),X(IBUF),
     &               X(LAI),X(IPAI),X(IPBJ),NOC,NVIR,NBF,X(LVEC),
     &               X(IERR),X(LC1),X(IAM),X(IX2),X(IPVT),X(IB),
     &               X(IODIIS),X(IENG),NROT,MAXIT)
C
      END IF
C
  450 CONTINUE
C
      CALL RETFM(NEEDF)
C
      CALL TSECND(TIM4)
      TIME = TIM4 - TIM3
      IF(MASWRK) WRITE(IW,9215) TIME
C
      IF(EXETYP.EQ.CHECK) GO TO 540
C
C     --- TRANSFORM COMPLETED CI-SINGLES DIFFERENCE DENSITY
C         TO THE AO BASIS AND WRITE TO DAF ---
C
      CALL CISAOD(X(LPIJ),X(LPAB),X(LAI),X(LVEC),X(LWRK1),X(LWRK2),
     &            X(LPMN),NBF,NBF2,NOC,NVIR,.TRUE.)
      CALL DAWRIT(IDAF,IODA,X(LPMN),NBF2,307,0)
C
C     --- ADD CI-SINGLES AO DIFFERENCE DENSITY TO HF DENSITY ---
C
      CALL DAREAD(IDAF,IODA,X(LWRK1),NBF2,16,0)
      CALL VADD(X(LPMN),1,X(LWRK1),1,X(LWRK2),1,NBF2)
      CALL DAWRIT(IDAF,IODA,X(LWRK1),NBF2,308,0)
      CALL DAWRIT(IDAF,IODA,X(LWRK2),NBF2,16,0)
C
C     --- FORM OCC-OCC, VIRT-VIRT, AND VIRT-OCC ENERGY WEIGHTED
C         DENSITY CONTRIBUTIONS ---
C
      CALL WOCPD(X(LWIJ),X(LPIJ),X(IENG),X(LS1),X(LPMN),X(LWRK1),
     &           X(LWRK2),X(LWRK3),X(NBUF),X(IBUF),X(LVEC),
     &           NOC,NBF,NBF2,NINTMX)
      CALL WVIRT(X(LWAB),X(LPAB),X(IENG),X(LS2),NOC,NVIR,NBF)
      CALL WVIROC(X(LAI),X(LWAI),X(IENG),X(LC2),NOC,NVIR,NBF)
C
C     --- TRANSFORM WIJ, WAB, AND WAI TO W(MUNU) I.E. MO -> AO ---
C
      CALL CISAOD(X(LWIJ),X(LWAB),X(LWAI),X(LVEC),X(LWRK1),X(LWRK2),
     &            X(LPMN),NBF,NBF2,NOC,NVIR,.TRUE.)
C
C     --- FORM TOTAL ENERGY WEIGHTED DENSITY MATRIX (HF+CIS)
C         AND WRITE TO DAF ---
C
      CALL WTOT(X(LVEC),X(IENG),X(LPMN),NOC,NBF,TWO)
      CALL DAWRIT(IDAF,IODA,X(LPMN),NBF2,309,0)
C
      CALL TSECND(TIM5)
      TIME = TIM5 - TIM4
      IF(MASWRK) WRITE(IW,9220) TIME
C
  540 CONTINUE
C
      CALL FLSHBF(IW)
C
      CALL RETFM(NEEDD)
      CALL RETFM(NEEDC)
      CALL RETFM(NEEDA)
C
C     --- DETERMINE CIS NATURAL ORBITALS ---
C
      IF(NPRINT.NE.-5) CALL CISNOS
C
      NT = NTSAVE
C
      RETURN
C
 9000 FORMAT(/10X,45('-')/
     &        13X,'ATOMIC ORBITAL BASIS CI-SINGLES ENERGY'/
     &        13X,'PROGRAM WRITTEN BY SIMON P. WEBB'/
     &        10X,45('-')/)
 9001 FORMAT(/10X,51('-')/
     &        10X,'ATOMIC ORBITAL BASIS CI-SINGLES ENERGY + PROPERTIES'/
     &        10X,'PROGRAM WRITTEN BY SIMON P. WEBB'/
     &        10X,51('-')/)
 9002 FORMAT(/10X,40('-')/
     &        10X,'ATOMIC ORBITAL BASIS CI-SINGLES GRADIENT'/
     &        14X,'PROGRAM WRITTEN BY SIMON P. WEBB'/
     &        10X,40('-')/)
 9005 FORMAT(1X,'# CORE ORBITALS      =',I5/
     &       1X,'# OCCUPIED ORBITALS  =',I5/
     &       1X,'# MOLECULAR ORBITALS =',I5/
     &       1X,'# BASIS FUNCTIONS    =',I5)
 9006 FORMAT(/1X,'NUMBER OF CIS SPIN-ADAPTED ANTISYMMETRIZED '
     &        ,'PRODUCTS (SAPS) IS',I8)
 9007 FORMAT(/1X,'NUMBER OF CIS DETERMINANTS IS',I8)
 9010 FORMAT(/1X,'NSTATE IS EQUAL TO THE NUMBER OF CIS CONFIGS'
     &       ,/1X,'SETTING MXVEC = NCFG'/)
 9015 FORMAT(/1X,'MXVEC CANNOT BE BIGGER THAN THE NUMBER OF CIS CONFIGS'
     &       ,/1X,'SETTING MXVEC = NCFG'/)
 9100 FORMAT(/5X,61('-')/
     &        7X,'SOLVE THE COUPLED-PERTURBED HARTREE-FOCK EQUATIONS'
     &           ,' IN THE',/,26X,'ATOMIC-ORBITAL BASIS'/
     &        5X,61('-')/)
 9200 FORMAT(/' -CIS- HAM. DIAGONAL ELEMENTS  TOOK ',F10.3,' SECONDS')
 9205 FORMAT(/' -CIS- ENERGY                  TOOK ',F10.3,' SECONDS')
 9210 FORMAT( ' -CIS- LAGRANGIAN              TOOK ',F10.3,' SECONDS')
 9215 FORMAT(/' -CIS- CPHF EQUATIONS          TOOK ',F10.3,' SECONDS')
 9220 FORMAT( ' -CIS- ENERGY WEIGHTED DENSITY TOOK ',F10.3,' SECONDS')
 9295 FORMAT(/1X,'APPROXIMATING CIS HAMILTONIAN DIAGONAL '
     &      ,'ELEMENTS USING ONLY ORBITAL ENERGIES')
 9297 FORMAT(/1X,'REGENERATING AO INTEGRAL LIST IN C1 POINT GROUP')
 9300 FORMAT(/1X,45('-')/' CALCULATING CIS HAMILTONIAN DIAGONAL '
     &      ,'ELEMENTS',/1X,45('-')/)
 9310 FORMAT(1X,'MIN MEMORY REQUIRED FOR FOCK-LIKE BUILDS =',I10,
     &       ' WORDS')
 9320 FORMAT(1X,'MEMORY REQ. FOR SINGLE BATCH BUILDS      =',I10,
     &       ' WORDS')
 9330 FORMAT(1X,'MEMORY AVAILABLE',25X,'=',I10,' WORDS')
 9340 FORMAT(/1X,'CARRYING OUT SINGLE BATCH CALCULATION')
 9350 FORMAT(/1X,'CARRYING OUT ',I5,' BATCH CALCULATION')
 9400 FORMAT(/1X,'MIN MEMORY REQ. FOR CIS ENERGY FOCK-LIKE BUILDS ='
     &       ,I10,' WORDS')
 9410 FORMAT(1X,'MEMORY REQ. FOR SINGLE BATCH BUILDS             ='
     &       ,I10,' WORDS')
 9420 FORMAT(1X,'MEMORY AVAILABLE',32X,'=',I10,' WORDS')
 9430 FORMAT(/1X,'SINGLE BATCH ENERGY CALCULATION WILL BE PERFORMED')
 9440 FORMAT(/1X,I5,' BATCH ENERGY CALCULATION WILL BE PERFORMED')
C
      END
C
C*MODULE CISGRD  *DECK CISAOE
      SUBROUTINE CISAOE(C,HFOE,BUF,IBUF,DIACIS,X,Y,A,B,
     &                  H1,H2,EIGV,VEC,SCR1,SCR2,ISCR,RNORM,
     &                  TA,TB,FA,FB,ESCF,NBF,NBF3,NOC,NVIR,
     &                  MXV,MXV2,NINTMX,NSTATE,ITST,DAVCVG,NDAVIT,
     &                  IA,EGY,SCLINT,NCFG,DETRUN,NCOR,NOCA,NGSVEC,
     &                  UNVGSS,NBATCH,NVBTCH,DGAPRX,IW,IP,MASWRK)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DOORTH,UNVGSS,DETRUN,DGAPRX,MASWRK
C
      DIMENSION C(NBF,NBF),HFOE(NBF),BUF(NINTMX),IBUF(NINTMX)
      DIMENSION DIACIS(NCFG),X(NBF,NBF),Y(NBF,NBF)
      DIMENSION A(NCFG,MXV),B(NCFG,MXV),H1(MXV2),H2(MXV2)
      DIMENSION EIGV(MXV),VEC(MXV*MXV),SCR1(8*MXV),SCR2(MXV)
      DIMENSION ISCR(MXV),RNORM(NSTATE),IA(*)
      DIMENSION TA(NBF,NBF,*),TB(NBF,NBF,*),FA(NBF*NBF,*),FB(NBF*NBF,*)
C
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00)
C
C     --- ROUTINE USING A DAVIDSON DIAGONALIZATION SCHEME TO FIND ---
C         EIGENVECTORS AND EIGENVALUES OF THE CI-SINGLES HAMILTONIAN
C         MATRIX IN THE ATOMIC ORBITAL BASIS
C
C         WRITTEN BY SIMON P. WEBB, JUNE 2003.
C
C     REQUIRED ON ENTRY:
C     -A-      INITIAL GUESS CIS VECTORS (IF INITIAL GUESS NOT PROVIDED
C              UNIT VECTORS ARE USED)
C     -C-      MO COEFFICENTS
C     -DIACIS- DIAGONAL ELEMENTS OF CIS HAMILTONIAN MATRIX
C
C     ON EXIT:
C     -A-     NSTATE CIS VECTORS
C     -EIGV-  NSTATE CIS ENERGIES
C
C     --- INITIALIZE SOME VARIABLES ---
C
      NGS2 = (NGSVEC*NGSVEC+NGSVEC)/2
      NROT = NOCA*NVIR
C
      M1 = NSTATE
      M2 = NSTATE
      MNXT = M2 + NSTATE
      IST = 1
C
C     --- GUESS UNIT VECTORS ORDERED ACCORDING TO CIS HAMILTONIAN
C         DIAGONAL ELEMENTS ---
C
      IF(UNVGSS) THEN
C
         IF(MASWRK) WRITE(IW,4000)
         CALL DCOPY(NCFG,DIACIS,1,B,1)
         CALL VCLR(A,1,NGSVEC*NCFG)
C
         IF(DETRUN) THEN
            COEF = ONE/SQRT(TWO)
            DO IIST=1,NGSVEC,2
               EE = 1.0D+08
               DO I=1,NROT
                  IF(B(I,1).LE.EE) THEN
                     EE = B(I,1)
                     ISAV1 = I
                  END IF
               END DO
               A(ISAV1,IIST) = COEF
               A(ISAV1+NROT,IIST) = -COEF
               A(ISAV1,IIST+1) = COEF
               A(ISAV1+NROT,IIST+1) = COEF
               B(ISAV1,1) = 1.0D+09
            END DO
         ELSE
            DO IIST=1,NGSVEC
               EE = 1.0D+08
               DO I=1,NROT
                  IF(B(I,1).LE.EE) THEN
                     EE = B(I,1)
                     ISAV = I
                  END IF
               END DO
               A(ISAV,IIST) = ONE
               B(ISAV,1) = 1.0D+09
            END DO
         END IF
C
      ELSE
         IF(MASWRK) WRITE(IW,4010)
      END IF
C
C     --- FIND INITIAL SET OF B = H*A VECTORS ---
C
      CALL HACIS(C,HFOE,BUF,IBUF,X,Y,A,B,TA,TB,FA,FB,
     &           NBF,NBF3,NOC,NVIR,NROT,NCFG,NINTMX,1,NGSVEC,
     &           SCLINT,DETRUN,NCOR,NOCA,NBATCH,NVBTCH)
C
C     --- FORM THE HAMILTONIAN IN THE VECTOR SUB-SPACE ---
C
      CALL VCLR(H1,1,MXV2)
C
      CALL HSUBSP(A,B,H1,NCFG,1,NGSVEC,EIGV,.FALSE.,IA)
      CALL DCOPY(NGS2,H1,1,H2,1)
C
C     --- DIAGONALIZE SMALL MATRIX ---
C
      IF(NGSVEC.EQ.NCFG) THEN
         IF(MASWRK) WRITE(IW,4400)
      END IF
      CALL GLDIAG(NGSVEC,NGSVEC,NGSVEC,H2,SCR1,EIGV,VEC,IGERR,ISCR)
      IF(IGERR.NE.0) CALL ABRT
C
C     --- FOR CASE WHERE FULL SPACE SPANNED BY NGSVEC WE ARE DONE ---
C
      IF(NGSVEC.EQ.NCFG) THEN
         CALL TRFRMB(A,NCFG,VEC,SCR2,NGSVEC,NGSVEC)
         DO IST=1,NSTATE
            IF(MASWRK) WRITE(IW,4410) IST,EIGV(IST)+ESCF
            IF(IST.EQ.ITST) EGY = EIGV(IST)+ESCF
         END DO
         IF(MASWRK) WRITE(IW,5045)
         RETURN
      END IF
C
C     --- TRANSFORM -A- AND -B- IN PREP FOR DAVIDSON ITERS ---
C
      IF(NGSVEC.GT.NSTATE) THEN
         CALL TRFRMB(A,NCFG,VEC,SCR2,NGSVEC,NSTATE)
         CALL TRFRMB(B,NCFG,VEC,SCR2,NGSVEC,NSTATE)
C
C        --- PUT FIRST NSTATE EIGV ELEMENTS ON DIAGONAL OF H1
C            AND PREPARE VEC ---
C
         CALL VCLR(H1,1,MXV2)
         CALL HSUBSP(A,B,H1,NCFG,1,NSTATE,EIGV,.TRUE.,IA)
C
         CALL VCLR(VEC,1,MXV*MXV)
         DO I = 1,NSTATE
            VEC((I-1)*NSTATE + I) = ONE
         END DO
      END IF
C
C         -------------------------------------------
C     --- START DAVIDSON ITERATIONS FOR TARGET STATES ---
C         -------------------------------------------
C
      IF(MASWRK) THEN
         WRITE(IW,4500)
         WRITE(IW,4510) NSTATE,NGSVEC,MXV,NDAVIT,DAVCVG
         WRITE(IW,4520)
         WRITE(IW,5000)
      END IF
C
      DO 100 ITGD=1,NDAVIT
C
C        --- IF MAXIMUM NUMBER OF A VECTORS REACHED THEN TRANSFORM
C            A AND B AND START WITH NEW INITIAL SET ---
C
         IF(MNXT.GT.MXV) THEN
            CALL TRFRMB(A,NCFG,VEC,SCR2,M2,NSTATE)
            CALL TRFRMB(B,NCFG,VEC,SCR2,M2,NSTATE)
            M1 = NSTATE
            M2 = NSTATE
C
C           --- PUT FIRST NSTATE EIGV ELEMENTS ON DIAGONAL OF H1
C               AND PREPARE VEC ---
C
            CALL VCLR(H1,1,MXV2)
            CALL HSUBSP(A,B,H1,NCFG,1,NSTATE,EIGV,.TRUE.,IA)
C
            CALL VCLR(VEC,1,MXV*MXV)
            DO I = 1,NSTATE
               VEC((I-1)*M2 + I) = ONE
            END DO
C
         END IF
C
C        --- FORM RESIDUE (GRADIENT) VECTORS AND AND THEIR NORMS ---
C
         ICOR = M1 + 1
         CALL RESVEC(A(1,ICOR),A,B,VEC,EIGV,NCFG,1,M2,RNORM,NSTATE)
C
         RNMRMS = DDOT(NSTATE,RNORM,1,RNORM,1)
         RNMRMS = SQRT(RNMRMS)
C
         DO IST=1,NSTATE
            IF(MASWRK) WRITE(IW,5020) IST,ITGD,EIGV(IST)+ESCF,
     &                                RNORM(IST)
         END DO
         IF(MASWRK) WRITE(IW,5045)
         IF(MASWRK) CALL FLSHBF(IW)
C
C        --- TEST NORM OF RESIDUE VECTOR FOR CONVERGENCE ---
C
         IF(RNMRMS.LT.DAVCVG .AND. ITGD.GT.1) THEN
C
            IF(MASWRK) WRITE(IW,5040) ITGD
C
            DO IST=1,NSTATE
               IF(MASWRK) WRITE(IW,5050) IST,EIGV(IST)+ESCF
               IF(IST.EQ.ITST) EGY = EIGV(IST)+ESCF
            END DO
            IF(MASWRK  .AND.  NSTATE.GT.1) WRITE(IW,5045)
C
            CALL TRFRMB(A,NCFG,VEC,SCR2,M2,NSTATE)
            CALL TRFRMB(B,NCFG,VEC,SCR2,M2,NSTATE)
C
            RETURN
         END IF
C
C       --- FORM FIRST ORDER CORRECTION FOR TARGET STATE EIGENVECTOR ---
C
         DO 200 IST=1,NSTATE
            ICOR = M1 + IST
            DO 210 II=1,NCFG
C
               VAL1 = EIGV(IST) - DIACIS(II)
               VAL2 = ABS(VAL1)
C
               IF(VAL2.LT.1.0D-12) THEN
                  VAL1 = 1.0D-05
               ELSE IF(VAL2.LT.1.0D-05) THEN
                  SIGN = VAL2/VAL1
                  VAL1 = SIGN*1.0D-05
               END IF
C
               A(II,ICOR) = A(II,ICOR)/VAL1
C
  210       CONTINUE
  200    CONTINUE
C
C        --- ORTHONORMALIZE TO MAKE NEW CORRECTION VECTOR ---
C
         DO 220 IST=1,NSTATE
C
            ICOR = M1 + IST
            DOORTH = .TRUE.
C
            DO 230 IORTH=1,5
               IF(.NOT.DOORTH) GO TO 230
               DOORTH = .FALSE.
               BN1 = DDOT(NCFG,A(1,ICOR),1,A(1,ICOR),1)
C
C              --- MAKE ORTHOGONAL TO ALL PREVIOUS VECTORS ---
C
               DO 250 M=1,ICOR-1
                  SUM = DDOT(NCFG,A(1,ICOR),1,A(1,M),1)
                  CALL DAXPY(NCFG,-SUM,A(1,M),1,A(1,ICOR),1)
  250          CONTINUE
C
C              --- CHECK WHETHER TO REPEAT ORTHOG. AND NORMALIZE ---
C
               BN2 = DDOT(NCFG,A(1,ICOR),1,A(1,ICOR),1)
C
               IF(BN2.LT.0.1D+00*BN1) DOORTH = .TRUE.
C
               SCALE = ONE/SQRT(BN2)
               CALL DSCAL(NCFG,SCALE,A(1,ICOR),1)
  230       CONTINUE
C
            IF(DOORTH) THEN
               IF(MASWRK)
     &         WRITE(IW,*)'ORTHOGANLIZATION PROBLEM IN CISAOE - ABORT'
               CALL ABRT
            END IF
C
  220    CONTINUE
C
C        --- ACT WITH H ON NEW VECTORS ---
C
         M2 = M2 + NSTATE
C
         CALL HACIS(C,HFOE,BUF,IBUF,X,Y,A,B,TA,TB,FA,FB,
     &              NBF,NBF3,NOC,NVIR,NROT,NCFG,NINTMX,M1+1,M2,
     &              SCLINT,DETRUN,NCOR,NOCA,NBATCH,NVBTCH)
C
C        --- FORM HAMILTONIAN IN NEW SUBSPACE BASIS ---
C
         CALL HSUBSP(A,B,H1,NCFG,M1+1,M2,EIGV,.FALSE.,IA)
         ICNT = (M2*M2+M2)/2
         CALL DCOPY(ICNT,H1,1,H2,1)
C
C        --- DIAGONALIZE SMALL MATRIX ---
C
         CALL GLDIAG(M2,M2,M2,H2,SCR1,EIGV,VEC,IGERR,ISCR)
         IF(IGERR.NE.0) CALL ABRT
C
         M1 = M2
         MNXT = M2 + NSTATE
C
  100 CONTINUE
C
      IF(MASWRK) THEN
         WRITE(IW,5100)
         IF(DGAPRX) WRITE(IW,5105)
         CALL CISVPU(IP,DETRUN,NSTATE,NCFG,EIGV,A)
      END IF
      CALL ABRT
C
      RETURN
C
 4000 FORMAT(/1X,'UNIT VECTOR GUESS AT CIS COEFFICIENTS ...')
 4010 FORMAT(/1X,'CIS GUESS COEFFICIENTS READ IN FROM INPUT',
     &        1X,'OR TAKEN FROM PREVIOUS STEP ...')
 4400 FORMAT(/1X,'USING IN MEMORY DIAGONALIZTION TO FIND CIS',
     &        1X,'EIGENVALUES AND EIGENVECTORS ...'/)
 4410 FORMAT(1X,'STATE',I5,' ENERGY=',F20.10)
 4500 FORMAT(/1X,65('-')/' USING DAVIDSON ALGORITHM TO FIND CIS ',
     &           'EIGENVALUES AND EIGENVECTORS'/,1X,65('-')/)
 4510 FORMAT(1X,'NUMBER OF STATES REQUESTED = ',I5/
     &       1X,'NUMBER OF GUESS VECTORS    = ',I5/
     &       1X,'MAX. NUMB. OF EXPAN. VECS. = ',I5/
     &       1X,'MAX. NUMB. OF ITERATIONS   = ',I5/
     &       1X,'CONVERGENCE CRITERION      = ',1P,E10.1)
 4520 FORMAT(/1X,'STARTING DAVIDSON ITERATIONS')
 5000 FORMAT(/1X,'STATE',6X,'ITERATION',6X,'ENERGY',9X,'RESIDUE NORM')
 5020 FORMAT(1X,I3,8X,I5,F20.10,F15.8)
 5040 FORMAT(/1X,'ALL STATES CONVERGED IN',1X,I4,1X,'ITERATIONS'/)
 5045 FORMAT(1X)
 5050 FORMAT(1X,'CONVERGED STATE',I5,' ENERGY=',F20.10)
 5100 FORMAT(1X,'ONE OR MORE STATES FAILED TO CONVERGE',/
     &       1X,'WRITING VECTORS TO .DAT AND ABORTING',/)
 5105 FORMAT(1X,'IF CONVERGING SLOWLY TRY SETTING DGAPRX=.FALSE.'/)
      END
C
C*MODULE CISGRD  *DECK HSUBSP
      SUBROUTINE HSUBSP(A,B,HMAT,NCFG,M1,M2,EIGV,DIAG,IA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DIAG
C
      DIMENSION A(NCFG,*),B(NCFG,*),HMAT(*),EIGV(*),IA(*)
C
      PARAMETER (ZERO=0.0D+00)
C
      IF(DIAG) THEN
         DO 10 N1 = M1,M2
            DO 20 N2 = 1,N1
               IDX = IA(N1) + N2
               IF(N2.EQ.N1) THEN
                  HMAT(IDX) = EIGV(N1)
               ELSE
                  HMAT(IDX) = ZERO
               END IF
   20       CONTINUE
   10    CONTINUE
C
      ELSE
C
         DO 50 N1 = M1,M2
            DO 60 N2 = 1,N1
C
               IDX = IA(N1) + N2
C
               DO 70 II = 1,NCFG
                  HMAT(IDX) = HMAT(IDX) + A(II,N2)*B(II,N1)
   70          CONTINUE
C
   60       CONTINUE
   50    CONTINUE
C
      END IF
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK RESVEC
      SUBROUTINE RESVEC(Q,A,B,VEC,EIGV,NCFG,M1,M2,RNORM,NSTATE)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION Q(NCFG,*),A(NCFG,*),B(NCFG,*),EIGV(*),VEC(M2,*)
      DIMENSION RNORM(NSTATE)
C
      PARAMETER (ZERO=0.0D+00)
C
      CALL VCLR(Q,1,NCFG*NSTATE)
C
      DO 5 IST=1,NSTATE
C
         DUM = ZERO
C
         DO 10 II = 1,NCFG
            DO 20 M = M1,M2
               Q(II,IST) = Q(II,IST) + VEC(M,IST)*(B(II,M)
     &                               - EIGV(IST)*A(II,M))
   20       CONTINUE
            DUM = DUM + Q(II,IST)*Q(II,IST)
   10    CONTINUE
C
         RNORM(IST) = SQRT(DUM)
C
    5 CONTINUE
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK TRFRMB
      SUBROUTINE TRFRMB(A,NCFG,VEC,WRK,M2,NST)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION A(NCFG,*),VEC(M2,*),WRK(*)
C
      PARAMETER (ZERO=0.0D+00)
C
      DO 10 II = 1,NCFG
C
         DO 20 JJ = 1,NST
            WRK(JJ) = ZERO
            DO 30 KK = 1,M2
               WRK(JJ) = WRK(JJ) + A(II,KK)*VEC(KK,JJ)
   30       CONTINUE
   20    CONTINUE
C
         DO 40 KI = 1,NST
            A(II,KI) = WRK(KI)
   40    CONTINUE
C
   10 CONTINUE
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK HACIS
      SUBROUTINE HACIS(C,HFOE,BUF,IBUF,X,Y,A,B,TA,TB,FA,FB,
     &           NBF,NBF3,NOC,NVIR,NROT,NCFG,NINTMX,NST1,NST2,
     &           SCAL,DETRUN,NCOR,NOCA,NBATCH,NVBTCH)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DETRUN
C
      DIMENSION C(NBF,NBF),HFOE(NBF),BUF(NINTMX),IBUF(NINTMX)
      DIMENSION X(NBF,NBF),Y(NBF,NBF),A(NCFG,*),B(NCFG,*)
      DIMENSION TA(NVBTCH*NBF3),TB(NVBTCH*NBF3)
      DIMENSION FA(NVBTCH*NBF3),FB(NVBTCH*NBF3)
C
C     --- VECTOR MULTIPLICATION OF TRIAL CI-SINGLES VECTOR -A-
C     BY CI-SINGLES HAMILTONIAN (B = H * A) - DO FOR EACH STATE ---
C
C     THIS IS DONE IN THE AO BASIS AVOIDING A 4-INDEX MO TRANSFORMATION
C     (FORESMAN ET AL. JPC(1992),VOL 96,PAGE 135 - SECTION 3)
C
      NV = NVBTCH
      NLOOP = NBATCH
      NVECS = NST2 - NST1 + 1
      NREM = 0
C
      IF(NVECS.LE.NVBTCH) THEN
         NV = NVECS
         NLOOP = 1
      END IF
C
C     --- FORM TRANSITION DENSITY -T-, CORRSPONDING FOCK-LIKE
C     MATRIX -F-, AND TRANSFORM -F- TO VIR-OCC MO BLOCK TO GIVE -B- ---
C
      DO 100 II=1,NLOOP
C
         IVEC = (II-1)*NV + 1
         IF(IVEC.GT.NVECS .OR. NREM.NE.0) GO TO 150
         IST = (II-1)*NV + NST1
         NREM = NVECS - IVEC + 1
         IF(NREM.LT.NV) THEN
            NV = NREM
         ELSE
            NREM = 0
         END IF
C
         CALL TFBAO(C,BUF,IBUF,X,Y,A(1,IST),B(1,IST),TA,TB,FA,FB,
     &              NCFG,NBF,NBF3,NOC,NVIR,NROT,NINTMX,SCAL,DETRUN,
     &              NCOR,NOCA,NCOR,NOCA,NV)
C
  100 CONTINUE
C
  150 CONTINUE
C
C     --- ADD HF ORBITAL ENERGY CONTRIBUTIONS ---
C
      DO 200 IST = NST1,NST2
C
         JJ = 0
         DO 210 I=NCOR+1,NOC
            DO 220 IDX=1,NVIR
C
               JJ = JJ + 1
               NA = IDX + NOC
               B(JJ,IST) = B(JJ,IST) + (HFOE(NA) - HFOE(I))
     &                               * A(JJ,IST)
C
               IF(DETRUN) THEN
                  B(JJ+NROT,IST) = B(JJ+NROT,IST) + (HFOE(NA) - HFOE(I))
     &                                            * A(JJ+NROT,IST)
               END IF
C
  220       CONTINUE
  210    CONTINUE
C
  200 CONTINUE
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK CISDGL
      SUBROUTINE CISDGL(DIACIS,C,X,FLM,PMN,XLMN,BUF,IBUF,HFOE,
     &                  NINTMX,NBF,NBF2,NOC,NVIR,SCAL,NCFG,
     &                  NCOR,NOCA,IA,NBATCH,NVBTCH,APPRX)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,APPRX
C
      DIMENSION DIACIS(NVIR,NOCA),C(NBF,NBF),X(NBF,NBF)
      DIMENSION FLM(NBF2),PMN(NVBTCH*NBF2),XLMN(NVBTCH*NBF2)
      DIMENSION BUF(NINTMX),IBUF(NINTMX),HFOE(NBF),IA(*)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (TWO=2.0D+00)
C
C     --- FORM THE DIAGONAL ELEMENTS OF THE CIS HAMILTONIAN MATRIX
C         IN THE AO BASIS ---
C
      CALL VCLR(DIACIS,1,NCFG)
C
      IF(APPRX) GO TO 200
C
      NV = NVBTCH
C
      IPCOUNT = ME - 1
C
C     --- LOOP OVER OCCUPIED ORBITAL BATCHES ---
C
      DO 10 IBTCH=1,NBATCH
C
         IF(IBTCH.EQ.NBATCH) NV = NOCA - ((NBATCH-1)*NVBTCH)
C
         DO 50 II=1,NV
C
C           --- FORM DENSITY MATRIX FOR ACTIVE OCCUPIED MO I ---
C
            I = (IBTCH-1)*NVBTCH + II + NCOR
C
            MUNU = 0
            DO 60 MU=1,NBF
               DUM = TWO*C(MU,I)
               DO 70 NU=1,MU
                  MUNU = MUNU + 1
                  IDX = (MUNU-1)*NV + II
                  PMN(IDX) = DUM*C(NU,I)
   70          CONTINUE
   60       CONTINUE
C
   50    CONTINUE
C
C        --- FORM FOCK-LIKE MATRIX FOR EACH OF ACTIVE MOS ---
C        DIRECT METHOD = RECOMPUTE 2E- AO INTEGRALS
C        STANDARD METHOD = PROCESS INTEGRALS FROM DISK
C
         CALL FLMSCL(PMN,XLMN,BUF,IBUF,NBF,NBF2,SCAL,NV)
C
         DO 150 II=1,NV
C
C           --- FIRST GATHER ELEMENTS OF THE FOCK-LIKE MATRIX
C
            CALL DCOPY(NBF2,XLMN(II),NV,FLM,1)
C
C           --- SYMMETRIZE EACH FOCK-LIKE MATRIX ---
C
            CALL SYMH(FLM,X,IA)
C
C           --- TRANSFORM BACK TO MO BASIS TO GIVE CONTRIBUTION TO
C               CIS HAMILTONIAN DIAGONAL ELEMENT ---
C
            JJ = (IBTCH-1)*NVBTCH + II
C
            DO 160 IDX=1,NVIR
               NA = IDX + NOC
C
               IF(GOPARR) THEN
                  IPCOUNT = IPCOUNT + 1
                  IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 160
               END IF
C
C              --- NOTE: OFF-DIAGONAL ELEMENTS OF -FLM- CALCULATED IN
C                        -FLMSCL- ARE DOUBLE THEIR TRUE VALUES ---
C
               LASI = 0
               DO 170 LAMD=1,NBF
                  DO 180 NSIG=1,LAMD
                     LASI = LASI + 1
C
                     DIACIS(IDX,JJ) = DIACIS(IDX,JJ) - C(LAMD,NA)*
     &                                C(NSIG,NA)*FLM(LASI)
C
  180             CONTINUE
  170          CONTINUE
  160       CONTINUE
C
  150    CONTINUE
C
   10 CONTINUE
C
C     --- SUM UP PARTIAL CONTRIBUTIONS IF PARALLEL ---
C
      IF(GOPARR) THEN
         CALL DDI_GSUMF(2316,DIACIS,NOCA*NVIR)
      END IF
C
C     --- ADD CONTRIBUTIONS FROM HF ORBITAL ENERGIES ---
C
  200 CONTINUE
C
      DO 210 II=1,NOCA
         I = II + NCOR
         DO 220 IDX=1,NVIR
C
            NA = IDX + NOC
            DIACIS(IDX,II) = DIACIS(IDX,II) + HFOE(NA) - HFOE(I)
C
  220    CONTINUE
  210 CONTINUE
C
      RETURN
      END
C
C*MODULE SCFLIB  *DECK FLMSCL
      SUBROUTINE FLMSCL(PMN,XLMN,BUF,IBUF,NBF,NBF2,SCAL,NV)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL PACK2E,SCHWRZ,DIRSCF,FDIFF,GOPARR,DSKWRK,MASWRK,TDSKWRK
C
      DIMENSION PMN(NV,NBF2),XLMN(NV,NBF2)
      DIMENSION BUF(NINTMX),IBUF(NINTMX)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH,
     &           MXAO=8192)
C
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     &                                                    NAUXSHL
      COMMON /FMCOM / XX(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     &                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     &                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     &                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (HALF=0.5D+00, TWO=2.0D+00)
C
      DATA RHF/8HRHF     /
C
C     --- FORM FOCK-LIKE MATRIX ---
C     DIRECT METHOD = RECOMPUTE 2E- AO INTEGRALS
C     STANDARD METHOD = PROCESS INTEGRALS FROM DISK
C
      NINT = 0
      NSCHWZ = 0
C
      IF(DIRSCF) THEN
C
C        --- HACK TO USE DFT SET UP TO SCALE EXCHANGE IN DIRFCK  ---
C
         SAV1 = DFTTYP(1)
         DFTTYP(1) = 1.0D+00
         SAV2 = DFTTYP(3)
         DFTTYP(3) = TWO*SCAL
C
         SCHWRZ=ISCHWZ.GT.0
         NSH2 = (NSHELL*NSHELL+NSHELL)/2
C
         CALL BASCHK(LMAX)
                       NANGM =  4
         IF(LMAX.EQ.2) NANGM =  6
         IF(LMAX.EQ.3) NANGM = 10
         IF(LMAX.EQ.4) NANGM = 15
         IF(LMAX.EQ.5) NANGM = 21
         IF(LMAX.EQ.6) NANGM = 28
         MAXG = NANGM**4
C
         CALL VALFM(LOADFM)
         LGHOND = LOADFM + 1
         LXINTS = LGHOND + MAXG
         LDSH   = LXINTS + NSH2
         LDDIJ  = LDSH   + NSH2
         LAST   = LDDIJ  + 49*MXG2
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
C
         IF(SCHWRZ) THEN
            DUMMY = 0.0D+00
            CALL AOSHLD(PMN,XX(LDSH),IA,NBF,NBF2,NSH2,NV)
            CALL DAREAD(IDAF,IODA,XX(LXINTS),NSH2,54,0)
         END IF
         CALL VCLR(XLMN,1,NV*NBF2)
         CALL TWOEI(RHF,DIRSCF,.FALSE.,.FALSE.,.FALSE.,
     &              INTTYP,SCHWRZ,NINT,NSCHWZ,NBF,NBF2,
     &              DUMMY,DUMMY,IDUMM,1,XX(LXINTS),NSH2,XX(LGHOND),
     &              MAXG,XX(LDDIJ),IA,PMN,XLMN,DUMMY,DUMMY,XX(LDSH),
     &              DUMMY,DUMMY,NV)
         CALL DSCAL(NBF2*NV,HALF,XLMN,1)
C
         CALL RETFM(NEED)
C
         DFTTYP(1) = SAV1
         DFTTYP(3) = SAV2
C
      ELSE
         TDSKWRK = DSKWRK
         DSKWRK  = .TRUE.
         CALL SEQREW(IJK)
         CALL HSTSCL(PMN,XLMN,BUF,IBUF,NINTMX,IA,NOPK,.FALSE.,SCAL,
     &               NBF2,NV)
         DSKWRK  = TDSKWRK
      END IF
C
C     --- SUM UP PARTIAL FOCK-LIKE MATRICES ---
C
      IF(GOPARR) THEN
         CALL DDI_GSUMF(2310,XLMN,NV*NBF2)
         CALL DDI_GSUMI(2311,NINT,1)
         CALL DDI_GSUMI(2312,NSCHWZ,1)
      END IF
C
      RETURN
      END
C
C*MODULE RHFUHF  *DECK HSTSCL
      SUBROUTINE HSTSCL(D,F,XX,IX,NINTMX,IA,NOPK,DUPAO,SCAL,NBF2,NV)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,DUPAO
C
      DIMENSION D(NV,NBF2),F(NV,NBF2),IA(*),XX(NINTMX),IX(*)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /ORDOPT/ NORDER(7)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
C
      PARAMETER (HALF=0.5D+00)
C
C     -- FORM A BATCH OF CLOSED SHELL FOCK-LIKE MATRICES ---
C     ONLY SYMMETRY UNIQUE 2E- INTEGRALS ARE PROCESSED, SO
C     THE RESULTING MATRIX MUST BE SYMMETRIZED BY -SYMH-
C
C     THE EXCHANGE INTEGRAL CONTRIBUTION IS SCALED ACCORDING ACCORDING
C     TO THE VALUE OF -SCAL-
C
C     INDICES IN AO INTEGRAL LABELS ARE IN STANDARD ORDER,
C               I.GE.J , K.GE.L , (IJ).GE.(KL)
C
C     IN CASE A DUPLICATED AO INTEGRAL LIST, AS OPPOSED TO A DISTRIBUTED
C     LIST, IS TO BE PRCESSED, WE MUST SKIP RECORDS IN ORDER TO MAKE THE
C     CPU WORK OF A FOCK BUILD BE PARALLELIZED.  DUPAO=.TRUE. IN CALL.
C
      CALL SEQREW(IS)
C
      CALL VCLR(F,1,NBF2*NV)
      I = 0
      J = 0
      K = 0
      L = 0
      NXX = 0
C
      IF(DUPAO) IPCOUNT = ME
C
C     --- INTEGRALS MUST BE IN STANDARD -J- LIST (NOPK=.TRUE.) ---
C
      IF(NOPK.NE.1) THEN
         WRITE(IW,*)'INTEGRALS MUST BE IN STANDARD -J- LIST'
         CALL ABRT
      END IF
C
  210 CONTINUE
C
C     DUPLICATED FILES MUST PREAD ONLY IF ITS OUR TURN TO DO A RECORD.
C     DISTRIBUTED FILES PROCESS EVERY RECORD IN THE LOCAL AO SUBFILE.
C
      IF(DUPAO) THEN
         IPCOUNT = IPCOUNT + 1
         IF (MOD(IPCOUNT,NPROC).NE.0) THEN
            READ (IS) NXX
            GO TO 280
         END IF
      END IF
C
      CALL PREAD(IS,XX,IX,NXX,NINTMX)
C
      IF(NXX .EQ. 0) GO TO 290
      NINT = IABS(NXX)
      IF(NINT .GT. NINTMX) CALL ABRT
C
      DO 260 M = 1,NINT
         NPACK = M
         IF (LABSIZ .EQ. 2) THEN
*I32        LABEL = IX( 2*NPACK - 1 )
*I32        IPACK = ISHFT( LABEL, -16 )
*I32        JPACK = IAND(  LABEL, 65535 )
*I32        LABEL = IX( 2*NPACK     )
*I32        KPACK = ISHFT( LABEL, -16 )
*I32        LPACK = IAND(  LABEL, 65535 )
*I64        LABEL = IX(NPACK)
*I64        IPACK = ISHFT( LABEL, -48 )
*I64        JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64        KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64        LPACK = IAND( LABEL, 65535 )
         ELSE IF (LABSIZ .EQ. 1) THEN
*I32        LABEL = IX(NPACK)
*I32        IPACK = ISHFT( LABEL, -24 )
*I32        JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32        KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I32        LPACK = IAND( LABEL, 255 )
*I64        IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64           LABEL = IX( NPACK/2 )
*I64           IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64           JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64           KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64           LPACK = IAND( LABEL, 255 )
*I64        ELSE
*I64           LABEL = IX( (NPACK/2)+1 )
*I64           IPACK = ISHFT( LABEL, -56 )
*I64           JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64           KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64           LPACK = IAND( ISHFT( LABEL, -32 ), 255 )
*I64        END IF
         END IF
C
         I = IPACK
         J = JPACK
         K = KPACK
         L = LPACK
C
         VAL = XX(M)
         NIJ = IA(I)+J
         NKL = IA(K)+L
C
C           USING SQUARE CANONICAL INTEGRAL FILE
C
         IF(NORDER(7).EQ.1) THEN
            IF(NKL .GT. NIJ) GO TO 260
            IF(I  .EQ.  J) VAL=VAL*HALF
            IF(K  .EQ.  L) VAL=VAL*HALF
            IF(NIJ.EQ.NKL) VAL=VAL*HALF
         END IF
C
         NIK = IA(I)+K
         NIL = IA(I)+L
         IF(J .LT. K) GO TO 220
         NJK = IA(J)+K
         NJL = IA(J)+L
         GO TO 240
C
  220    NJK = IA(K)+J
         IF(J .LT. L) GO TO 230
         NJL = IA(J)+L
         GO TO 240
C
  230    NJL = IA(L)+J
C
  240    CONTINUE
C
C        --- SCALE ACCORDING TO HAMTYP AND MULT ---
C
         VAL2 = (VAL+VAL)
         VALS = VAL*SCAL
C
C        --- FORM CONTRIBUTIONS FOR BATCH OF MOS ---
C
         DO IMO=1,NV
            F(IMO,NIJ) = F(IMO,NIJ) + VAL2*D(IMO,NKL)
            F(IMO,NKL) = F(IMO,NKL) + VAL2*D(IMO,NIJ)
C
            F(IMO,NIK) = F(IMO,NIK) - VALS*D(IMO,NJL)
            F(IMO,NIL) = F(IMO,NIL) - VALS*D(IMO,NJK)
            F(IMO,NJK) = F(IMO,NJK) - VALS*D(IMO,NIL)
            F(IMO,NJL) = F(IMO,NJL) - VALS*D(IMO,NIK)
         END DO
C
  260 CONTINUE
C
C     DONE PROCESSING THIS ENTIRE INTEGRAL RECORD, GO BACK FOR MORE
C
  280 CONTINUE
      IF(NXX.GT.0) GO TO 210
C
C     DONE CONSTRUCTING MATRIX, NOTE: OFF DIAGONAL ELEMENTS ARE DOUBLE
C
  290 CONTINUE
C
      CALL SEQREW(IS)
C
      RETURN
      END
C
C*MODULE RHFUHF  *DECK HSTNSY
      SUBROUTINE HSTNSY(DA,DB,FA,FB,XX,IX,NINTMX,IA,NOPK,DUPAO,NBF,
     &                  SCAL,DETRUN,NV)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DUPAO,DETRUN
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION DA(NV,NBF,NBF),DB(NV,NBF,NBF)
      DIMENSION FA(NV,NBF,NBF),FB(NV,NBF,NBF)
      DIMENSION XX(NINTMX),IX(*),IA(*)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /ORDOPT/ NORDER(7)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
C
      PARAMETER (HALF=0.5D+00)
C
C     -- FORM BATCH OF NON-SYMMETRIC CLOSED SHELL FOCK-LIKE MATRICES --
C
C     INDICES IN AO INTEGRAL LABELS ARE IN STANDARD ORDER,
C               I.GE.J , K.GE.L , (IJ).GE.(KL)
C     THIS ROUTINE ALSO HAS AN OPTION TO PROCESS A -P- INTEGRAL LIST,
C     AND THERE IS SPECIAL CODE FOR FORMING A MATRIX WITH A DESIRED
C     ADMIXTURE OF THE EXCHANGE TERM (FOR PURPOSES OF DFT MATRICES).
C
C     IN CASE A DUPLICATED AO INTEGRAL LIST, AS OPPOSED TO A DISTRIBUTED
C     LIST, IS TO BE PRCESSED, WE MUST SKIP RECORDS IN ORDER TO MAKE THE
C     CPU WORK OF A FOCK BUILD BE PARALLELIZED.  DUPAO=.TRUE. IN CALL.
C
      CALL SEQREW(IS)
C
      I = 0
      J = 0
      K = 0
      L = 0
      NXX = 0
C
      IF(DUPAO) IPCOUNT = ME
C
C     ----- INTEGRALS MUST BE IN STANDARD -J- LIST (NOPK=.TRUE.) -----
C
      IF(NOPK.NE.1) THEN
         WRITE(IW,*)'INTEGRALS MUST BE IN STANDARD -J- LIST'
         CALL ABRT
      END IF
C
  210 CONTINUE
C
C     DUPLICATED FILES MUST PREAD ONLY IF ITS OUR TURN TO DO A RECORD.
C     DISTRIBUTED FILES PROCESS EVERY RECORD IN THE LOCAL AO SUBFILE.
C
      IF(DUPAO) THEN
         IPCOUNT = IPCOUNT + 1
         IF (MOD(IPCOUNT,NPROC).NE.0) THEN
            READ (IS) NXX
            GO TO 280
         END IF
      END IF
C
      CALL PREAD(IS,XX,IX,NXX,NINTMX)
C
      IF(NXX .EQ. 0) GO TO 290
      NINT = IABS(NXX)
      IF(NINT .GT. NINTMX) CALL ABRT
C
      DO 260 M = 1,NINT
         NPACK = M
         IF (LABSIZ .EQ. 2) THEN
*I32        LABEL = IX( 2*NPACK - 1 )
*I32        IPACK = ISHFT( LABEL, -16 )
*I32        JPACK = IAND(  LABEL, 65535 )
*I32        LABEL = IX( 2*NPACK     )
*I32        KPACK = ISHFT( LABEL, -16 )
*I32        LPACK = IAND(  LABEL, 65535 )
*I64        LABEL = IX(NPACK)
*I64        IPACK = ISHFT( LABEL, -48 )
*I64        JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64        KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64        LPACK = IAND( LABEL, 65535 )
         ELSE IF (LABSIZ .EQ. 1) THEN
*I32        LABEL = IX(NPACK)
*I32        IPACK = ISHFT( LABEL, -24 )
*I32        JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32        KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I32        LPACK = IAND( LABEL, 255 )
*I64        IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64           LABEL = IX( NPACK/2 )
*I64           IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64           JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64           KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64           LPACK = IAND( LABEL, 255 )
*I64        ELSE
*I64           LABEL = IX( (NPACK/2)+1 )
*I64           IPACK = ISHFT( LABEL, -56 )
*I64           JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64           KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64           LPACK = IAND( ISHFT( LABEL, -32 ), 255 )
*I64        END IF
         END IF
C
         I = IPACK
         J = JPACK
         K = KPACK
         L = LPACK
C
         VAL = XX(M)
         NIJ = IA(I)+J
         NKL = IA(K)+L
C
C        --- MAY BE USING SQUARE CANONICAL INTEGRAL FILE ---
C
         IF(NORDER(7).EQ.1) THEN
            IF(NKL .GT. NIJ) GO TO 260
            IF(I  .EQ.  J) VAL=VAL*HALF
            IF(K  .EQ.  L) VAL=VAL*HALF
            IF(NIJ.EQ.NKL) VAL=VAL*HALF
         END IF
C
C        --- SCALE ACCORDING TO HAMTYP AND MULT ---
C
         VAL2 = VAL*SCAL
C
C     --- LOOP OVER BATCH OF MOS TO CONTRIBUTE TO FOCK-LIKE MATRIX ---
C
         DO IMO=1,NV
            DUMLK = VAL2*(DA(IMO,K,L)+DA(IMO,L,K))
            DUMIJ = VAL2*(DA(IMO,I,J)+DA(IMO,J,I))
            FA(IMO,I,J) = FA(IMO,I,J) + DUMLK
            FA(IMO,J,I) = FA(IMO,J,I) + DUMLK
            FA(IMO,K,L) = FA(IMO,K,L) + DUMIJ
            FA(IMO,L,K) = FA(IMO,L,K) + DUMIJ
C
            FA(IMO,I,K) = FA(IMO,I,K) - VAL*DA(IMO,J,L)
            FA(IMO,K,I) = FA(IMO,K,I) - VAL*DA(IMO,L,J)
            FA(IMO,I,L) = FA(IMO,I,L) - VAL*DA(IMO,J,K)
            FA(IMO,L,I) = FA(IMO,L,I) - VAL*DA(IMO,K,J)
            FA(IMO,J,K) = FA(IMO,J,K) - VAL*DA(IMO,I,L)
            FA(IMO,K,J) = FA(IMO,K,J) - VAL*DA(IMO,L,I)
            FA(IMO,J,L) = FA(IMO,J,L) - VAL*DA(IMO,I,K)
            FA(IMO,L,J) = FA(IMO,L,J) - VAL*DA(IMO,K,I)
         END DO
C
C        --- FOR DETERMINANT BASED RUN BOTH ALPHA AND BETA
C        MATRICES REQUIRED ---
C
         IF(DETRUN) THEN
            DO IMO=1,NV
               DUMLK = VAL2*(DA(IMO,K,L)+DA(IMO,L,K))
               DUMIJ = VAL2*(DA(IMO,I,J)+DA(IMO,J,I))
               FB(IMO,I,J) = FB(IMO,I,J) + DUMLK
               FB(IMO,J,I) = FB(IMO,J,I) + DUMLK
               FB(IMO,K,L) = FB(IMO,K,L) + DUMIJ
               FB(IMO,L,K) = FB(IMO,L,K) + DUMIJ
C
               DUMLK = VAL2*(DB(IMO,K,L)+DB(IMO,L,K))
               DUMIJ = VAL2*(DB(IMO,I,J)+DB(IMO,J,I))
               FB(IMO,I,J) = FB(IMO,I,J) + DUMLK
               FB(IMO,J,I) = FB(IMO,J,I) + DUMLK
               FB(IMO,K,L) = FB(IMO,K,L) + DUMIJ
               FB(IMO,L,K) = FB(IMO,L,K) + DUMIJ
C
               FA(IMO,I,J) = FA(IMO,I,J) + DUMLK
               FA(IMO,J,I) = FA(IMO,J,I) + DUMLK
               FA(IMO,K,L) = FA(IMO,K,L) + DUMIJ
               FA(IMO,L,K) = FA(IMO,L,K) + DUMIJ
C
               FB(IMO,I,K) = FB(IMO,I,K) - VAL*DB(IMO,J,L)
               FB(IMO,K,I) = FB(IMO,K,I) - VAL*DB(IMO,L,J)
               FB(IMO,I,L) = FB(IMO,I,L) - VAL*DB(IMO,J,K)
               FB(IMO,L,I) = FB(IMO,L,I) - VAL*DB(IMO,K,J)
               FB(IMO,J,K) = FB(IMO,J,K) - VAL*DB(IMO,I,L)
               FB(IMO,K,J) = FB(IMO,K,J) - VAL*DB(IMO,L,I)
               FB(IMO,J,L) = FB(IMO,J,L) - VAL*DB(IMO,I,K)
               FB(IMO,L,J) = FB(IMO,L,J) - VAL*DB(IMO,K,I)
            END DO
         END IF
C
  260 CONTINUE
C
C     DONE PROCESSING THIS ENTIRE INTEGRAL RECORD, GO BACK FOR MORE
C
  280 CONTINUE
      IF(NXX.GT.0) GO TO 210
  290 CONTINUE
C
      CALL SEQREW(IS)
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK CISAOL
      SUBROUTINE CISAOL(XLAI,PIJ,PAB,PMN,C,X,Y,XLMN,BUF,IBUF,
     *                  C1,C2,NOC,NVIR,NBF,NINTMX)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION XLAI(NVIR,NOC),PIJ(NOC,NOC),PAB(NVIR,NVIR),PMN(*)
      DIMENSION C(NBF,*),X(NBF,NBF),Y(NBF,NBF),XLMN(*)
      DIMENSION BUF(NINTMX),IBUF(NINTMX)
      DIMENSION C1(NVIR,NOC),C2(NVIR,NOC)
C
      PARAMETER (ONE=1.0D+00)
C
C     --- CONSTRUCTION OF THE CI-SINGLES LAGRANGIAN ---
C
C     INPUT:
C        C   : MO COEFFICIENTS
C        PIJ : OCC-OCC CIS DIFFERENCE DENSITY (MO BASIS)
C        PAB : VIRT-VIRT CIS DIFFERENCE DENSITY (MO BASIS)
C        C1  : INTERMEDIATE ARRAY FORMED IN SUBROUTINE -CISITM-
C        C2  : INTERMEDIATE ARRAY FORMED IN SUBROUTINE -CISITM-
C
C     OUTPUT:
C        XLAI : CIS LAGRANGIAN MATRIX (MO BASIS)
C
C     --- TRANSFORM -PIJ- AND -PAB- TO AO BASIS AS -PMN- ---
C
      CALL MRARTR(PIJ,NOC,NOC,NOC,C,NBF,NBF,X,NBF)
      CALL MRARBR(C,NBF,NBF,NOC,X,NBF,NBF,Y,NBF)
C
      CALL MRARTR(PAB,NVIR,NVIR,NVIR,C(1,NOC+1),NBF,NBF,X,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIR,ONE,C(1,NOC+1),NBF,X,NBF,
     &           ONE,Y,NBF)
C
      MUNU = 0
      DO 10 MU=1,NBF
         DO 20 NU=1,MU
            MUNU = MUNU + 1
            PMN(MUNU) = Y(MU,NU) + Y(NU,MU)
   20    CONTINUE
   10 CONTINUE
C
C     --- THIRD AND FOURTH LAGRANGIAN TERMS, IN AO BASIS ---
C         DIRECT METHOD = RECOMPUTE 2E- AO INTEGRALS
C         STANDARD METHOD = PROCESS INTEGRALS FROM DISK
C
      CALL FLMAT(PMN,Y,XLMN,BUF,IBUF,NBF)
C
C     --- TRANSFORM THESE TWO TERMS TO MO BASIS, AND ADD ALL
C         CONTRIBUTIONS TOGETHER ---
C
      CALL VCLR(XLAI,1,NOC*NVIR)
      CALL MRARBR(Y,NBF,NBF,NBF,C,NBF,NOC,X,NBF)
      CALL DGEMM('T','N',NVIR,NOC,NBF,ONE,C(1,NOC+1),NBF,X,NBF,
     &           ONE,XLAI,NVIR)
C
      DO 30 MI=1,NOC
         DO 40 MA=1,NVIR
            XLAI(MA,MI) = -XLAI(MA,MI) - C1(MA,MI) + C2(MA,MI)
   40    CONTINUE
   30 CONTINUE
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK CISITM
      SUBROUTINE CISITM(C,F,X,FOCC,FVIR,A,B,C1,C2,S1,S2,PIJ,PTEMP,
     &                  HFOE,NOC,NVIR,NBF,NCOR,NOCA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION C(NBF,*),F(NBF,NBF),X(*)
      DIMENSION FOCC(NOC,NOCA),FVIR(NVIR,NVIR)
      DIMENSION A(NVIR,NOCA),B(NVIR,NOC)
      DIMENSION C1(NVIR,NOC),C2(NVIR,NOC),S1(NOC,NOC),S2(NVIR,NVIR)
      DIMENSION PIJ(NOC,NOC),PTEMP(NOCA,NOCA),HFOE(NBF)
C
      IF(NCOR.GT.0) THEN
         CALL VCLR(C1,1,NCOR*NVIR)
         CALL VCLR(S1,1,NCOR*NOC)
      END IF
C
C     --- TRANSFORM FOCK-LIKE MATRIX F (FORMED USING THE TRANSITION
C         DENSITY) TO OCC-OCC AND VIR-VIR BLOCKS IN MO BASIS ---
C
      CALL MRARBR(F,NBF,NBF,NBF,C(1,NCOR+1),NBF,NOCA,X,NBF)
      CALL MRTRBR(C,NBF,NBF,NOC,X,NBF,NOCA,FOCC,NOC)
C
      CALL MRARBR(F,NBF,NBF,NBF,C(1,NOC+1),NBF,NVIR,X,NBF)
      CALL MRTRBR(C(1,NOC+1),NBF,NBF,NVIR,X,NBF,NVIR,FVIR,NVIR)
C
C     --- FORM C1 AND C2 ---
C
      CALL MRTRBR(FVIR,NVIR,NVIR,NVIR,A,NVIR,NOCA,C1(1,NCOR+1),NVIR)
      CALL MRARTR(A,NVIR,NVIR,NOCA,FOCC,NOC,NOC,C2,NVIR)
C
C     --- FORM S1 AND S2 ---
C
      CALL MRTRBR(B,NVIR,NVIR,NOC,A,NVIR,NOCA,S1(1,NCOR+1),NOC)
      CALL MRARTR(A,NVIR,NVIR,NOCA,B(1,NCOR+1),NVIR,NVIR,S2,NVIR)
C
C     --- FROZEN CORE-ACTIVE DENSITY MATRIX ELEMENTS ---
C
      IF(NCOR.GT.0) THEN
         DO 10 J=NCOR+1,NOC
            DO 20 I=1,NOC
               IF(I.LE.NCOR) THEN
                  EIJ = HFOE(J) - HFOE(I)
                  PIJ(I,J) = S1(I,J)/EIJ
                  PIJ(J,I) = PIJ(I,J)
               ELSE
                  PIJ(I,J) = PTEMP(I-NCOR,J-NCOR)
               END IF
   20       CONTINUE
   10    CONTINUE
      END IF
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK CISAOD
      SUBROUTINE CISAOD(PIJ,PAB,PAI,C,X,Y,PMN,NBF,NBF2,NOC,NVIR,
     &                  RESDEN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL RESDEN
C
      PARAMETER (ONE=1.0D+00)
C
      DIMENSION PIJ(NOC,NOC),PAB(NVIR,NVIR),PAI(NVIR,NOC)
      DIMENSION C(NBF,*),X(NBF,*),Y(NBF,*),PMN(*)
C
C     --- TRANSFORM VIR-OCC BLOCK TO AO BASIS AND SYMMETRIZE ---
C
      IF(RESDEN) THEN
C
         CALL MRARTR(PAI,NVIR,NVIR,NOC,C,NBF,NBF,X,NBF)
         CALL MRARBR(C(1,NOC+1),NBF,NBF,NVIR,X,NBF,NBF,Y,NBF)
C
         MUNU = 0
         DO 15 MU=1,NBF
            DO 10 NU=1,MU
               MUNU = MUNU + 1
               PMN(MUNU) = Y(MU,NU) + Y(NU,MU)
   10       CONTINUE
   15    CONTINUE
C
      ELSE
C
         CALL VCLR(PMN,1,NBF2)
C
      END IF
C
C     --- TRANSFORM OCC-OCC AND VIR-VIR BLOCKS TO AO BASIS ---
C
      CALL MRARTR(PIJ,NOC,NOC,NOC,C,NBF,NBF,X,NBF)
      CALL MRARBR(C,NBF,NBF,NOC,X,NBF,NBF,Y,NBF)
C
      CALL MRARTR(PAB,NVIR,NVIR,NVIR,C(1,NOC+1),NBF,NBF,X,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIR,ONE,C(1,NOC+1),NBF,X,NBF,ONE,
     &           Y,NBF)
C
C     --- COMPLETE DENSITY MATRIX IN AO BASIS (LOWER TRIANGLE) ---
C
      MUNU = 0
      DO 25 MU=1,NBF
         DO 20 NU=1,MU
            MUNU = MUNU + 1
            PMN(MUNU) = PMN(MUNU) + Y(MU,NU)
   20    CONTINUE
   25 CONTINUE
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK TFBAO
      SUBROUTINE TFBAO(C,BUF,IBUF,X,Y,A,B,TA,TB,FA,FB,NCFG,
     &                 NBF,NBF3,NOC,NVIR,NROT,NINTMX,SCAL,
     &                 DETRUN,NCOR,NOCA,NC,NO,NV)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DETRUN
C
      DIMENSION C(NBF,NBF),BUF(NINTMX),IBUF(NINTMX)
      DIMENSION X(NBF,NBF),Y(NBF,NBF),A(NCFG,*),B(NCFG,*)
      DIMENSION TA(NV*NBF3),TB(NV*NBF3),FA(NV*NBF3),FB(NV*NBF3)
C
C     --- TRANSFORM CIS VECTOR -A- TO AO BASIS GIVING
C     TRANSITION DENSITIY -T- ---
C
      DO 50 II=1,NV
C
         CALL MRARTR(A(1,II),NVIR,NVIR,NOCA,C(1,NCOR+1),NBF,NBF,X,NBF)
         CALL MRARBR(C(1,NOC+1),NBF,NBF,NVIR,X,NBF,NBF,Y,NBF)
         CALL DCOPY(NBF3,Y,1,TA(II),NV)
C
         IF(DETRUN) THEN
            CALL MRARTR(A(NROT+1,II),NVIR,NVIR,NOCA,C(1,NCOR+1),
     &                  NBF,NBF,X,NBF)
            CALL MRARBR(C(1,NOC+1),NBF,NBF,NVIR,X,NBF,NBF,Y,NBF)
            CALL DCOPY(NBF3,Y,1,TB(II),NV)
         END IF
C
   50 CONTINUE
C
C     --- FORM FOCK-LIKE MATRIX -F- USING TRANSITION DENSITY -T- ---
C
C     DIRECT METHOD = RECOMPUTE 2E- AO INTEGRALS
C     STANDARD METHOD = PROCESS INTEGRALS FROM DISK
C
      CALL VCLR(FA,1,NBF*NBF*NV)
      IF(DETRUN) CALL VCLR(FB,1,NBF*NBF*NV)
C
      CALL FLMNSY(TA,TB,FA,FB,BUF,IBUF,NBF,SCAL,DETRUN,NV)
C
C     --- TRANSFORM FOCK-LIKE MATRIX -F- TO MO BASIS
C     (VIRTUAL-OCCUPIED BLOCK ONLY) THIS IS -B- ---
C
      DO 100 II=1,NV
C
         CALL DCOPY(NBF3,FA(II),NV,Y,1)
         CALL MRARBR(Y,NBF,NBF,NBF,C(1,NC+1),NBF,NO,X,NBF)
         CALL MRTRBR(C(1,NOC+1),NBF,NBF,NVIR,X,NBF,NO,B(1,II),NVIR)
C
         IF(DETRUN) THEN
            CALL DCOPY(NBF3,FB(II),NV,Y,1)
            CALL MRARBR(Y,NBF,NBF,NBF,C(1,NCOR+1),NBF,NO,X,NBF)
            CALL MRTRBR(C(1,NOC+1),NBF,NBF,NVIR,X,NBF,NO,
     &                  B(NROT+1,II),NVIR)
         END IF
C
  100 CONTINUE
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK WVIRT
      SUBROUTINE WVIRT(WAB,PAB,E,S2,NOC,NVIR,NBF)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (PT5=0.5D+00)
C
      DIMENSION WAB(NVIR,NVIR),PAB(NVIR,NVIR),E(NBF),S2(NVIR,NVIR)
C
C     --- FORM VIRT-VIRT BLOCK OF THE ENERGY-WEIGHTED DENSITY ---
C
C     INPUT:
C     PAB : VIRT-VIRT DIFFERENCE DENSITY
C     E   : ORBITAL ENERGIES
C     S2  : INTERMEDIATE ARRAY FORMED IN SUBROUTINE -CISITM-
C
C     OUTPUT:
C     WAB : VIRT-VIRT ENERGY WEIGHTED DENSITY
C
      DO 50 IA=1,NVIR
         DO 60 IB=1,NVIR
            WAB(IB,IA) = -PAB(IB,IA)*(E(IB+NOC)+E(IA+NOC))*PT5
     &                   -(S2(IB,IA)+S2(IA,IB))*PT5
   60    CONTINUE
   50 CONTINUE
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK WOCPD
      SUBROUTINE WOCPD(WIJ,PIJ,E,S1,PMN,XLMN,X,Y,BUF,IBUF,
     &                 C,NOC,NBF,NBF2,NINTMX)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00, PT5=0.5D+00)
C
      DIMENSION WIJ(NOC,NOC),PIJ(NOC,NOC),E(*),S1(NOC,NOC)
      DIMENSION C(NBF,NBF),XLMN(NBF,NBF),PMN(NBF2)
      DIMENSION X(NBF,NBF),Y(NBF,NBF),BUF(NINTMX),IBUF(NINTMX)
C
C     --- FORM OCC-OCC BLOCK OF THE ENERGY-WEIGHTED DENSITY ---
C
C     INPUT:
C     PIJ : OCC-OCC CIS DIFFERENCE DENSITY IN MO BASIS
C     E   : ORBITAL ENERGIES
C     PMN : COMPLETE CIS DIFFERENCE DENSITY IN THE AO BASIS
C     S1  : INTERMEDIATE ARRAY FORMED IN SUBROUTINE -CISITM-
C
C     OUTPUT:
C     WIJ : OCC-OCC ENERGY WEIGHTED DENSITY IN MO BASIS
C
C     --- ORBITAL ENERGY CONTRIBUTIONS TO -WIJ- ---
C
      DO 25 J=1,NOC
         DO 20 I=1,NOC
            WIJ(I,J) = -PIJ(I,J)*(E(I)+E(J))*PT5
   20    CONTINUE
   25 CONTINUE
C
C     --- INTEGRAL TERM OF THE ENERGY WEIGHTED DENSITY IN AO BASIS ---
C
      CALL FLMAT(PMN,Y,XLMN,BUF,IBUF,NBF)
C
C     --- TRANSFORM THE FOCK-LIKE MATRIX TO MO BASIS AND ADD ---
C
      CALL MRARBR(Y,NBF,NBF,NBF,C,NBF,NOC,X,NBF)
      CALL DGEMM('T','N',NOC,NOC,NBF,-TWO,C,NBF,X,NBF,ONE,WIJ,NOC)
C
C     --- COMPLETE THE ENERGY WEIGHTED DENSITY IN THE MO BASIS ---
C
      DO 50 J=1,NOC
         DO 60 I=1,NOC
            WIJ(I,J) = WIJ(I,J) - (S1(I,J)+S1(J,I))*PT5
   60    CONTINUE
   50 CONTINUE
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK WVIROC
      SUBROUTINE WVIROC(PAI,WAI,E,C2,NOC,NVIR,NBF)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION PAI(NVIR,NOC),WAI(NVIR,NOC),E(NBF),C2(NVIR,NOC)
C
C     --- FORM VIRT-OCC BLOCK OF THE ENERGY-WEIGHTED DENSITY ---
C
C     INPUT:
C     PAI : VIRT-OCC CIS DIFFERENCE DENSITY IN MO BASIS
C     E   : ORBITAL ENERGIES
C     C2  : INTERMEDIATE ARRAY FORMED IN SUBROUTINE -CISITM-
C
C     OUTPUT:
C     WAI : VIRT-OCC ENERGY WEIGHTED DENSITY IN MO BASIS
C
      DO 10 I=1,NOC
         DO 20 NA=1,NVIR
            WAI(NA,I) = -E(I)*PAI(NA,I) - C2(NA,I)
   20    CONTINUE
   10 CONTINUE
C
      RETURN
      END
C
C*MODULE MP2GRD  *DECK JKDCIS
      SUBROUTINE JKDCIS(PHF,PCIS,TCIS,NBF3,NBF2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
C
      DIMENSION PHF(NBF2),PCIS(NBF2),TCIS(NBF2)
C
C     --- READ ONE-PARTICLE SCF DENSITY, CIS DENSITY DIFFERENCE,
C         CIS TRANSITION DENSITY MATRICES, AND MOS ---
C
      CALL DAREAD(IDAF,IODA,PHF ,NBF2,308,0)
      CALL DAREAD(IDAF,IODA,PCIS,NBF2,307,0)
      CALL DAREAD(IDAF,IODA,TCIS,NBF3,301,0)
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK WTOT
      SUBROUTINE WTOT(C,E,W2,NOC,NBF,OCC)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00)
C
      DIMENSION W2(*),C(NBF,*),E(*)
C
C     ---  W2 IS THE HF+CIS ENERGY-WEIGHTED DENSITY  ---
C
      IJ = 0
      DO 25 I=1,NBF
         DO 20 J=1,I
            IJ = IJ + 1
            DUM = ZERO
            DO 30 K=1,NOC
               DUM = DUM - E(K)*C(I,K)*C(J,K)
   30       CONTINUE
            W2(IJ) = W2(IJ) + DUM*OCC
   20    CONTINUE
   25 CONTINUE
C
      RETURN
      END
C
C*MODULE SCFLIB  *DECK FLMNSY
      SUBROUTINE FLMNSY(DA,DB,FA,FB,BUF,IBUF,NBF,SCAL,DETRUN,NV)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,TDSKWRK
      LOGICAL PACK2E,SCHWRZ,DIRSCF,FDIFF,DETRUN,DIRCIS,DET
C
      DIMENSION DA(NV*NBF*NBF),DB(NV*NBF*NBF)
      DIMENSION FA(NV*NBF*NBF),FB(NV*NBF*NBF)
      DIMENSION BUF(NINTMX),IBUF(NINTMX)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH,
     &           MXAO=8192)
C
      COMMON /CISSCL/ SCLINT,DET
      COMMON /FMCOM / XX(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     &                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     &                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     &                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA RHF/8HRHF     /
C
C     --- FORM SQUARE NON-SYMMETRIC FOCK-LIKE MATRIX ---
C         DIRECT METHOD = RECOMPUTE 2E- AO INTEGRALS
C         STANDARD METHOD = PROCESS INTEGRALS FROM DISK
C
      SCLINT = SCAL
      DET = DETRUN
      NINT = 0
      NSCHWZ = 0
      NBF2 = (NBF*NBF+NBF)/2
      NBF3 = NBF*NBF
C
      IF(DIRSCF) THEN
         SCHWRZ=ISCHWZ.GT.0
         NSH2 = (NSHELL*NSHELL+NSHELL)/2
C
         CALL BASCHK(LMAX)
                       NANGM =  4
         IF(LMAX.EQ.2) NANGM =  6
         IF(LMAX.EQ.3) NANGM = 10
         IF(LMAX.EQ.4) NANGM = 15
         IF(LMAX.EQ.5) NANGM = 21
         IF(LMAX.EQ.6) NANGM = 28
         MAXG = NANGM**4
C
         CALL VALFM(LOADFM)
         LGHOND = LOADFM + 1
         LXINTS = LGHOND + MAXG
         LDSH   = LXINTS + NSH2
         LDDIJ  = LDSH   + NSH2
         LAST   = LDDIJ  + 49*MXG2
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
C
         IF(SCHWRZ) THEN
            DUMMY = 0.0D+00
            CALL SHLSQD(RHF,DA,DB,DUMMY,XX(LDSH),IA,NBF,NBF2,
     &                  NSH2,1,NV,DETRUN)
            CALL DAREAD(IDAF,IODA,XX(LXINTS),NSH2,54,0)
         END IF
         DIRCIS = .TRUE.
         CALL TWOEI(RHF,.FALSE.,.FALSE.,.FALSE.,DIRCIS,
     &              INTTYP,SCHWRZ,NINT,NSCHWZ,NBF,NBF2,
     &              DUMMY,DUMMY,IDUMM,1,XX(LXINTS),NSH2,XX(LGHOND),
     &              MAXG,XX(LDDIJ),IA,DA,FA,DB,FB,XX(LDSH),
     &              DUMMY,DUMMY,NV)
C
         CALL RETFM(NEED)
      ELSE
         TDSKWRK = DSKWRK
         DSKWRK  = .TRUE.
         CALL SEQREW(IJK)
         CALL HSTNSY(DA,DB,FA,FB,BUF,IBUF,NINTMX,IA,NOPK,.FALSE.,NBF,
     &               SCAL,DETRUN,NV)
         DSKWRK  = TDSKWRK
      END IF
C
C     --- SUM UP PARTIAL FOCK-LIKE MATRICES ---
C
      IF(GOPARR) THEN
         CALL DDI_GSUMI(2311,NINT,1)
         CALL DDI_GSUMI(2312,NSCHWZ,1)
         CALL DDI_GSUMF(2313,FA,NV*NBF3)
         IF(DETRUN) CALL DDI_GSUMF(2314,FB,NV*NBF3)
      END IF
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK CISSPN
      SUBROUTINE CISSPN(A,NSPN,NROT,NCFG,NSTATE,DETRUN,MULT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DETRUN
C
      PARAMETER (ZERO=0.0D+00)
C
      DIMENSION A(NCFG,NSTATE),NSPN(NSTATE)
C
C     --- WORK OUT SPIN OF EACH STATE ---
C
      IF(DETRUN) THEN
C
         DO IST=1,NSTATE
            DUM1 = ZERO
            DO I=1,NROT
               DUM2 = ABS(A(I,IST))
               IF(DUM2.GT.DUM1) THEN
                  DUM1 = DUM2
                  IROT = I
               END IF
            END DO
            SIGN = A(IROT,IST)/A(IROT+NROT,IST)
            IF(SIGN.GT.ZERO) THEN
               NSPN(IST) = 0
            ELSE
               NSPN(IST) = 1
            END IF
         END DO
C
      ELSE
         DO IST=1,NSTATE
            IF(MULT.EQ.1) THEN
               NSPN(IST) = 0
            ELSE
               NSPN(IST) = 1
            END IF
         END DO
      END IF
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK CISCMB
      SUBROUTINE CISCMB(A,NSPN,NROT,NCFG,NSTATE,ISTATE,SCLINT,CISSNG)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL CISSNG
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
C
      DIMENSION A(NCFG,NSTATE),NSPN(NSTATE)
C
C     --- TAKE APPROPRIATE LINEAR COMBINATIONS OF CIS COEFFICENTS
C         FOR SPIN OF TARGET STATE ISTATE ---
C
      COEF = ONE/SQRT(TWO)
C
      IF(NSPN(ISTATE).EQ.0) THEN
C
C        --- SINGLET STATE ---
C
         SCLINT = TWO
         CISSNG = .TRUE.
         DO I=1,NROT
            A(I,ISTATE) = COEF*(A(I,ISTATE) + A(I+NROT,ISTATE))
         END DO
C
      ELSE IF(NSPN(ISTATE).EQ.1) THEN
C
C        --- TRIPLET STATE ---
C
         SCLINT = ZERO
         CISSNG = .FALSE.
         DO I=1,NROT
            A(I,ISTATE) = COEF*(A(I,ISTATE) - A(I+NROT,ISTATE))
         END DO
C
      ELSE
         WRITE(6,*)'PROBLEM IN CISCMB - ABORTING'
         CALL ABRT
      END IF
C
      RETURN
      END
C
C*MODULE CISGRD  *DECK CISOUT
      SUBROUTINE CISOUT(A,CISE,NSPN,ISVNAM,ESCF,NROT,NCFG,NSTATE,NCOR,
     &                  NOCA,NOC,NVIR,PRTTOL,DETRUN,IW)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL DETRUN,CHECK
C
      DIMENSION A(NCFG,NSTATE),CISE(NSTATE),NSPN(NSTATE),ISVNAM(NSTATE)
C
      COMMON /SOSYM/  EULANG(4,48),GAM(48,48),IRMON(MXAO)
      COMMON /SYMMUL/ NIJREP(14,14),IJREP(2,14,14,14)
      COMMON /SYMREP/ IRPNAM(14),IPA(14),LAMBDA(14),LAMBD0(14),
     &                IADDR1(14),IADDR2(14),IADDR3(14)
C
      PARAMETER (ONE=1.0D+00, ONED7=1.0D+07)
C
C     CONVERSION FACTORS FROM NATIONAL INSTITUTE OF STANDARDS AND TECH.
C
      PARAMETER (HTOEV=27.2113845D+00, HTOKLM=627.50956D+00)
      PARAMETER (HTOWN=2.194746313705D+05)
C
C     --- OUTPUT CI-SINGLES ENERGIES, SPINS, AND CI COEFFICIENTS ---
C
      IF(DETRUN) THEN
         WRITE(IW,7190)
         WRITE(IW,6000)
         WRITE(IW,7191)
      ELSE
         WRITE(IW,7190)
         WRITE(IW,6050)
         WRITE(IW,7191)
      END IF
C
      WRITE(IW,7000) PRTTOL
      WRITE(IW,7005) ESCF
C
      DO 10 IST=1,NSTATE
C
         NAMSAV = 1
         CHECK = .FALSE.
         ICFG = 0
         DO IFRM=1,NOCA
            DO ITO=1,NVIR
               ICFG = ICFG + 1
               IF(ABS(A(ICFG,IST)).GE.PRTTOL) THEN
                  MO1 = IFRM + NCOR
                  MO2 = ITO + NOC
                  IRP = IRMON(MO1)
                  JRP = IRMON(MO2)
                  IINAM = IJREP(2,1,IRP,JRP)
                  IF(CHECK .AND. IINAM.NE.NAMSAV) THEN
                     WRITE(IW,7008)
                     IINAM = 1
                     GO TO 50
                  END IF
                  NAMSAV = IINAM
                  CHECK = .TRUE.
               END IF
            END DO
         END DO
C
   50    CONTINUE
         ISVNAM(IST) = IINAM
C
         ESTATE = CISE(IST)+ESCF
         SPIN = NSPN(IST)
         WRITE(IW,7010) IST,ESTATE,SPIN,IRPNAM(IINAM)
C
C        --- ALPHA AND BETA DETERMINANTS ---
C
         IF(DETRUN) THEN
C
            WRITE(IW,7050)
            ICFG = 0
            DO IFRM=1,NOCA
               DO ITO=1,NVIR
                  ICFG = ICFG + 1
                  IF(ABS(A(ICFG,IST)).GE.PRTTOL) THEN
                     MO1 = IFRM + NCOR
                     MO2 = ITO + NOC
                     WRITE(IW,7070) MO1,MO2,A(ICFG,IST)
                  END IF
               END DO
            END DO
C
            ICFG = 0
            DO IFRM=1,NOCA
               DO ITO=1,NVIR
                  ICFG = ICFG + 1
                  IF(ABS(A(ICFG+NROT,IST)).GE.PRTTOL) THEN
                     MO1 = IFRM + NCOR
                     MO2 = ITO + NOC
                     WRITE(IW,7072) MO1,MO2,A(ICFG+NROT,IST)
                  END IF
               END DO
            END DO
            WRITE(IW,7195)
C
C        --- SPIN-ADAPTED ANTISYMMETRIZED PRODUCTS ---
C
         ELSE
C
            WRITE(IW,7150)
            WRITE(IW,7200)
            ICFG = 0
            DO IFRM=1,NOCA
               DO ITO=1,NVIR
                  ICFG = ICFG + 1
                  IF(ABS(A(ICFG,IST)).GE.PRTTOL) THEN
                     MO1 = IFRM + NCOR
                     MO2 = ITO + NOC
                     WRITE(IW,7160) MO1,MO2,A(ICFG,IST)
                  END IF
               END DO
            END DO
            WRITE(IW,7200)
C
         END IF
   10 CONTINUE
C
C     --- OUTPUT TABLE OF EXCITATION ENERGIES ---
C
      WRITE(IW,8000)
      DO IST=1,NSTATE
         IINAM = ISVNAM(IST)
         MULT = 2*NSPN(IST) + 1
         WRITE(IW,8010) MULT,IRPNAM(IINAM),CISE(IST),
     &   CISE(IST)*HTOEV,CISE(IST)*HTOKLM,CISE(IST)*HTOWN,
     &   (ONE/(CISE(IST)*HTOWN))*ONED7
      END DO
C
      RETURN
 6000 FORMAT(8X,'RESULTS FROM DETERMINANT BASED ATOMIC ORBITAL'
     &       ,' CI-SINGLES',/,28X,'ENERGY CALCULATION')
 6050 FORMAT(8X,'RESULTS FROM SPIN-ADAPTED ANTISYMMETRIZED PRODUCT',
     &       ' (SAPS)',/,8X,'BASED ATOMIC ORBITAL CI-SINGLES ENERGY',
     &       ' CALCULATION')
 7000 FORMAT(/1X,'PRINTING CIS COEFFICIENTS LARGER THAN',F10.6)
 7005 FORMAT(/1X,'RHF REFERENCE ENERGY  = ',F20.10)
 7008 FORMAT(/1X,'WARNING: THIS STATE HAS BROKEN SYMMETRY, CHECK MOS')
 7010 FORMAT(//1X,'EXCITED STATE',I4,'  ENERGY= ',F20.10,'  S =',F4.1,
     &           '  SPACE SYM = ',A4)
 7050 FORMAT(/7X,53(1H-),/19X,'SINGLE EXCITATION',14X,
     &        'COEFFICENT',/,12X,'FROM SPIN ORB',5X,'TO SPIN ORB',/,
     &        7X,53(1H-))
 7070 FORMAT(7X,'ALPHA',3X,I5,12X,I5,10X,F12.8)
 7072 FORMAT(7X,'BETA ',3X,I5,12X,I5,10X,F12.8)
 7150 FORMAT(/10X,46(1H-),/10X,'SINGLE EXCITATION',15X,'SAP COEFFICENT',
     &       /,10X,'FROM MO',5X,'TO MO')
 7160 FORMAT(9X,I5,6X,I5,17X,F12.8)
 7190 FORMAT(/8X,56(1H-))
 7191 FORMAT(8X,56(1H-))
 7195 FORMAT(7X,53(1H-))
 7200 FORMAT(10X,46(1H-))
 8000 FORMAT(/2X,69(1H-),/,20X,'CI-SINGLES EXCITATION ENERGIES',
     &       /,2X,'STATE',7X,'HARTREE',8X,'EV',6X,'KCAL/MOL',7X,
     &       'CM-1',9X,'NM'/,2X,69(1H-))
 8010 FORMAT(3X,I1,A4,2X,F12.10,4X,F7.4,4X,F8.4,4X,F10.2,4X,F7.2)
      END
C
C*MODULE SCFLIB  *DECK SHLSQD
      SUBROUTINE SHLSQD(SCFTYP,DA,DB,DGVB,DSH,IA,L1,L2,NSH2,NFO,NV,
     &                  DETRUN)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL ALPHA,BETA,DETRUN
C
      DIMENSION DA(NV,L1,L1),DB(NV,L1,L1)
      DIMENSION DGVB(NV,L2,NFO),DSH(NSH2),IA(L1)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
C
      PARAMETER (ZERO=0.0D+00)
C
      DATA RHF,UHF,ROHF,GVB/8HRHF     ,8HUHF     ,8HROHF    ,8HGVB     /
C
C     --- COMPRESS DENSITY OVER AOS TO DENSITY OVER SHELLS
C     FOR THE CASE OF MULTIPLE SQUARE DENSITY MATRICES ---
C
      IJSH=0
      ALPHA = SCFTYP.EQ.RHF  .OR.  SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
      BETA  = DETRUN         .OR.  SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
      DO 240 ISH=1,NSHELL
         MINI = KLOC(ISH)
         MAXI = MINI + KMAX(ISH) - KMIN(ISH)
         DO 230 JSH=1,ISH
            MINJ = KLOC(JSH)
            MAXJ = MINJ + KMAX(JSH) - KMIN(JSH)
            IJSH = IJSH+1
            DMAX = ZERO
            DO 130 I=MINI,MAXI
               IF(ISH.EQ.JSH) MAXJ=I
               DO 120 J=MINJ,MAXJ
                  IJ = IA(I) + J
                  DO IV=1,NV
                     IF(ALPHA) THEN
                        IF(ABS(DA(IV,J,I)).GT.DMAX)
     &                     DMAX = ABS(DA(IV,J,I))
                     END IF
                     IF(BETA) THEN
                        IF(ABS(DB(IV,J,I)).GT.DMAX)
     &                     DMAX = ABS(DB(IV,J,I))
                     END IF
                     IF(SCFTYP.EQ.GVB) THEN
                        DO 110 IHAM=1,NFO
                           IF(ABS(DGVB(IV,IJ,IHAM)).GT.DMAX)
     *                                  DMAX=ABS(DGVB(IV,IJ,IHAM))
  110                   CONTINUE
                     END IF
                  END DO
  120          CONTINUE
  130       CONTINUE
            DSH(IJSH) = DMAX
  230    CONTINUE
  240 CONTINUE
      RETURN
      END
C
C*MODULE SCFLIB  *DECK DRFCIS
      SUBROUTINE DRFCIS(DA,FA,DB,FB,GHONDO,NBF,NINT,NV)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL IANDJ,KANDL,SAME,OUT,DETRUN
C
      DIMENSION GHONDO(*)
      DIMENSION DA(NV,NBF,NBF),DB(NV,NBF,NBF)
      DIMENSION FA(NV,NBF,NBF),FB(NV,NBF,NBF)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      COMMON /CISSCL/ SCAL,DETRUN
      COMMON /ERIOUT/ ISH,JSH,KSH,LSH,LSTRI,LSTRJ,LSTRK,LSTRL
      COMMON /MISC  / IANDJ,KANDL,SAME
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJX,IJ,KL,IJKL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
C
      DATA HALF /0.5D+00/
C
C-NEXT STATEMENT IS FOR VARIOUS IBM XLF 3.X AND 5.X COMPILERS-
C
      SAVE IJN,KLN
C
C     --- FORM A BATCH OF SQUARE NON-SYMMETRIC FOCK-LIKE MATRICES
C     DIRECTLY FROM INTEGRALS ---
C
      CUTINT = CUTOFF
C
      SAME  = ISH .EQ. KSH .AND. JSH .EQ. LSH
      IANDJ = ISH .EQ. JSH
      KANDL = KSH .EQ. LSH
C
      MINI = KMIN(ISH)
      MINJ = KMIN(JSH)
      MINK = KMIN(KSH)
      MINL = KMIN(LSH)
      MAXI = KMAX(ISH)
      MAXJ = KMAX(JSH)
      MAXK = KMAX(KSH)
      MAXL = KMAX(LSH)
      LOCI = KLOC(ISH)-MINI
      LOCJ = KLOC(JSH)-MINJ
      LOCK = KLOC(KSH)-MINK
      LOCL = KLOC(LSH)-MINL
C
      IJN = 0
      JMAX = MAXJ
      DO 360 I = MINI,MAXI
         I_ERIC = (I-MINI)*LSTRI + 1
         IF (IANDJ) JMAX = I
         DO 340 J = MINJ,JMAX
            IJ_ERIC = (J-MINJ)*LSTRJ + I_ERIC
            IJN = IJN+1
            LMAX = MAXL
            KLN = 0
            DO 320 K =  MINK,MAXK
               IJK_ERIC = (K-MINK)*LSTRK + IJ_ERIC
               IF (KANDL) LMAX = K
               DO 300 L = MINL,LMAX
                  KLN = KLN+1
                  IF(SAME .AND. KLN.GT.IJN) GO TO 340
C
                  IJKL_ERIC = (L-MINL)*LSTRL + IJK_ERIC
                  VAL = GHONDO( IJKL_ERIC )
                  IF(ABS(VAL).LT.CUTINT) GO TO 300
                  NINT = NINT + 1
C
                  II = LOCI+I
                  JJ = LOCJ+J
                  KK = LOCK+K
                  LL = LOCL+L
                  IF (II .GE. JJ) GO TO 100
                  N = II
                  II = JJ
                  JJ = N
  100             IF (KK .GE. LL) GO TO 120
                  N = KK
                  KK = LL
                  LL = N
  120             IF (II-KK) 140,160,180
  140             N = II
                  II = KK
                  KK = N
                  N = JJ
                  JJ = LL
                  LL = N
                  GO TO 180
  160             IF (JJ .LT. LL) GO TO 140
  180             CONTINUE
C
                  IF(II.EQ.JJ) VAL = VAL*HALF
                  IF(KK.EQ.LL) VAL = VAL*HALF
                  IF(II.EQ.KK  .AND.  JJ.EQ.LL) VAL = VAL*HALF
C
                  VAL2 = VAL*SCAL
C
C                 --- LOOP OVER BATCH OF MOS AND MAKE CONTRIB. TO
C                 FOCK-LIKE MATRIX ---
C
                  DO IMO=1,NV
                     DUMLK = VAL2*(DA(IMO,KK,LL)+DA(IMO,LL,KK))
                     DUMIJ = VAL2*(DA(IMO,II,JJ)+DA(IMO,JJ,II))
                     FA(IMO,II,JJ) = FA(IMO,II,JJ) + DUMLK
                     FA(IMO,JJ,II) = FA(IMO,JJ,II) + DUMLK
                     FA(IMO,KK,LL) = FA(IMO,KK,LL) + DUMIJ
                     FA(IMO,LL,KK) = FA(IMO,LL,KK) + DUMIJ
C
                     FA(IMO,II,KK) = FA(IMO,II,KK) - VAL*DA(IMO,JJ,LL)
                     FA(IMO,KK,II) = FA(IMO,KK,II) - VAL*DA(IMO,LL,JJ)
                     FA(IMO,II,LL) = FA(IMO,II,LL) - VAL*DA(IMO,JJ,KK)
                     FA(IMO,LL,II) = FA(IMO,LL,II) - VAL*DA(IMO,KK,JJ)
                     FA(IMO,JJ,KK) = FA(IMO,JJ,KK) - VAL*DA(IMO,II,LL)
                     FA(IMO,KK,JJ) = FA(IMO,KK,JJ) - VAL*DA(IMO,LL,II)
                     FA(IMO,JJ,LL) = FA(IMO,JJ,LL) - VAL*DA(IMO,II,KK)
                     FA(IMO,LL,JJ) = FA(IMO,LL,JJ) - VAL*DA(IMO,KK,II)
                  END DO
C
C                 --- FOR DETERMINANT BASED RUN BOTH ALPHA AND BETA
C                     MATRICES ARE REQUIRED ---
C
                  IF(DETRUN) THEN
                  DO IMO=1,NV
                     DUMLK = VAL2*(DA(IMO,KK,LL)+DA(IMO,LL,KK))
                     DUMIJ = VAL2*(DA(IMO,II,JJ)+DA(IMO,JJ,II))
                     FB(IMO,II,JJ) = FB(IMO,II,JJ) + DUMLK
                     FB(IMO,JJ,II) = FB(IMO,JJ,II) + DUMLK
                     FB(IMO,KK,LL) = FB(IMO,KK,LL) + DUMIJ
                     FB(IMO,LL,KK) = FB(IMO,LL,KK) + DUMIJ
C
                     DUMLK = VAL2*(DB(IMO,KK,LL)+DB(IMO,LL,KK))
                     DUMIJ = VAL2*(DB(IMO,II,JJ)+DB(IMO,JJ,II))
                     FB(IMO,II,JJ) = FB(IMO,II,JJ) + DUMLK
                     FB(IMO,JJ,II) = FB(IMO,JJ,II) + DUMLK
                     FB(IMO,KK,LL) = FB(IMO,KK,LL) + DUMIJ
                     FB(IMO,LL,KK) = FB(IMO,LL,KK) + DUMIJ
C
                     FA(IMO,II,JJ) = FA(IMO,II,JJ) + DUMLK
                     FA(IMO,JJ,II) = FA(IMO,JJ,II) + DUMLK
                     FA(IMO,KK,LL) = FA(IMO,KK,LL) + DUMIJ
                     FA(IMO,LL,KK) = FA(IMO,LL,KK) + DUMIJ
C
                     FB(IMO,II,KK) = FB(IMO,II,KK) - VAL*DB(IMO,JJ,LL)
                     FB(IMO,KK,II) = FB(IMO,KK,II) - VAL*DB(IMO,LL,JJ)
                     FB(IMO,II,LL) = FB(IMO,II,LL) - VAL*DB(IMO,JJ,KK)
                     FB(IMO,LL,II) = FB(IMO,LL,II) - VAL*DB(IMO,KK,JJ)
                     FB(IMO,JJ,KK) = FB(IMO,JJ,KK) - VAL*DB(IMO,II,LL)
                     FB(IMO,KK,JJ) = FB(IMO,KK,JJ) - VAL*DB(IMO,LL,II)
                     FB(IMO,JJ,LL) = FB(IMO,JJ,LL) - VAL*DB(IMO,II,KK)
                     FB(IMO,LL,JJ) = FB(IMO,LL,JJ) - VAL*DB(IMO,KK,II)
                  END DO
                  END IF
C
  300          CONTINUE
  320       CONTINUE
  340    CONTINUE
  360 CONTINUE
C
      RETURN
      END
C*MODULE CISGRD  *DECK CISNOS
      SUBROUTINE CISNOS
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      LOGICAL MNMEDG,MNMEOP,UNVGSS,DGAPRX,RDCISV
C
      PARAMETER (MXATM=2000)
C
      COMMON /CISEN / ECIS
      COMMON /CISPAR/ HAMTYP,DIAGZN,DAVCVG,PRTTOL,NSTATE,ISTATE,MULT,
     &                MXV,NDAVIT,ICISPR,NACORE,NBCORE,NOA,NOB,NORBOC,
     &                NBF,NGSVEC,MNMEDG,MNMEOP,ICLOBBR,UNVGSS,DGAPRX,
     &                RDCISV
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,MNUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
C     --- ALLOCATE MEMORY TO GENERATE CIS NATURAL ORBITALS ---
C
      L0 = NQMT
      L1 = NBF
      L2 = (NBF*NBF+NBF)/2
      L3 = NBF*NBF
      CALL VALFM(LOADFM)
      LEIG   = LOADFM + 1
      LIWRK  = LEIG   + L1
      LSCR   = LIWRK  + L1
      LDA    = LSCR   + 8*L1
      LDB    = LDA    + L2
      LVEC   = LDB    + L2
      LSPHLF = LVEC   + L3
      LSMHLF = LSPHLF + L3
      LAST   = LSMHLF + L3
      NEEDNO = LAST - LOADFM - 1
      CALL GETFM(NEEDNO)
      IF(EXETYP.EQ.CHECK) GO TO 800
C
      CALL UHFNOS(X(LEIG),X(LIWRK),X(LSCR),X(LDA),X(LDB),X(LVEC),
     *            X(LSPHLF),X(LSMHLF),L0,L1,L2,1)
C
      OCCSUM = 0.0D+00
      VIRSUM = 0.0D+00
      DO 100 I=1,NOA
         OCCSUM = OCCSUM + X(LEIG-1+I)
  100 CONTINUE
      DO 110 I=NOA+1,L0
         VIRSUM = VIRSUM + X(LEIG-1+I)
  110 CONTINUE
C
      IF(MASWRK) THEN
         WRITE(IW,9110) ISTATE
         WRITE(IW,9120) (X(LEIG-1+I),I=1,L0)
         WRITE(IW,9130) OCCSUM,VIRSUM
         WRITE(IW,9140)
         CALL PREV(X(LVEC),X(LEIG),L0,L1,L1)
         WRITE(IW,9180)
         WRITE(IP,8000) ECIS,ISTATE
         WRITE(IP,8010) ' $VEC'
         CALL PUSQL(X(LVEC),L0,L1,L1)
         WRITE(IP,8010) ' $END'
         CALL FLSHBF(IP)
      END IF
C
  800 CONTINUE
      CALL RETFM(NEEDNO)
      RETURN
C
 8000 FORMAT('CIS NATURAL ORBITALS, E(CIS)=',F20.10,2X,'EXCITED STATE '
     &       ,I5)
 8010 FORMAT(A5)
 9110 FORMAT(/1X,'CIS NATURAL ORBITAL OCCUPATION NUMBERS FOR EXCITED',
     *           ' STATE',I5,' ARE')
 9120 FORMAT(1X,10F7.4)
 9130 FORMAT(1X,'THERE ARE',F11.4,' ELECTRONS IN PRINCIPAL CIS',
     &          ' NATURAL ORBITALS.'/
     &       1X,'THERE ARE',F11.4,' ELECTRONS IN SECONDARY CIS',
     &          ' NATURAL ORBITALS.')
 9140 FORMAT(/10X,20(1H-)/10X,'CIS NATURAL ORBITALS'/
     *        10X,20(1H-))
 9180 FORMAT(/1X,'THE CIS NATURAL ORBITALS HAVE BEEN PUNCHED.')
      END
C*MODULE CISGRD  *DECK CISVPU
      SUBROUTINE CISVPU(IP,DETRUN,NSTATE,NCFG,ESTATE,CIVEC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION CIVEC(NCFG,NSTATE),ESTATE(NSTATE)
      LOGICAL DETRUN
      CHARACTER*4 BASIS
C
C        PUNCH CI VECTORS IN $CISVEC GROUP
C        THIS ROUTINE TO BE CALLED ONLY BY MASTER NODE
C
      WRITE(IP,FMT='(8H $CISVEC)')
                 BASIS='SAPS'
      IF(DETRUN) BASIS='DETS'
      WRITE(IP,9010) NCFG,BASIS,NSTATE
      DO IST=1,NSTATE
         WRITE(IP,9020) IST,ESTATE(IST)
         WRITE(IP,9030) (CIVEC(ICFG,IST),ICFG=1,NCFG)
      ENDDO
      WRITE(IP,FMT='(5H $END)')
C
      CALL FLSHBF(IP)
      RETURN
C
 9010 FORMAT(I10,1X,A4,',',I4,' STATES')
 9020 FORMAT('STATE',I4,' ENERGY=',F20.10)
 9030 FORMAT(1P,5E14.6)
      END
C*MODULE CISGRD  *DECK CISVRD
      SUBROUTINE CISVRD(DETRUN,NCFG,NSTATE)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DETRUN,GOPARR,DSKWRK,MASWRK,TDSKWRK
      CHARACTER*4 BASIS,BASFIL
C
      COMMON /CIFILS/ NFT11,NFT12,NFT13,NFT14,NFT15,NFT16,IDAF20,NEMEMX
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     --- READ IN CIS VECTORS FROM INPUT FILE AND WRITE TO DISK ---
C
      CALL VALFM(LOADFM)
      LVEC  = 1     + LOADFM
      LAST  = LVEC  + NCFG*NSTATE
      NEED  = LAST  - LOADFM
      CALL GETFM(NEED)
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $CISVEC',IEOF)
      IF (IEOF.EQ.1) THEN
         IF(MASWRK) WRITE(IW,*) '$CISVEC NOT FOUND'
         CALL ABRT
      END IF
C
C        READ HEADER VALUES, AND CHECK THEM AGAINST THIS RUN
C
                 BASIS='SAPS'
      IF(DETRUN) BASIS='DETS'
      IF(MASWRK) READ(IR,9010) NCFFIL,BASFIL,NSTFIL
      CALL DDI_BCAST(1027,'I',NCFFIL,1,MASTER)
      CALL DDI_BCAST(1028,'I',NSTFIL,1,MASTER)
      IF(NCFG.NE.NCFFIL  .OR.  NSTATE.GT.NSTFIL) THEN
         IF(MASWRK) WRITE(IW,9040) NCFG,BASIS,NCFFIL,BASFIL,
     *                             NSTATE,NSTFIL
         CALL ABRT
      END IF
C
C        READ THE STATES, WITHOUT ANY CHECKING, SKIPPING ENERGY LINES
C
      IF(MASWRK) THEN
         WRITE(IW,9000) NSTATE,NCFG
         LPOINT = LVEC-1
         DO I=1,NSTATE
            READ(IR,*)
            READ(IR,9030) (X(LPOINT+ICFG),ICFG=1,NCFG)
            LPOINT = LPOINT+NCFG
         ENDDO
      END IF
C
C     --- WRITE CIS VECTORS TO DISK FILE 12 ---
C     SINCE FILE 12 IS MAINTAINED ONLY ON THE MASTER NODE,
C     WE DO NOT NEED TO BROADCAST THE INPUT CI VECTORS NOW,
C     SO JUST WRITE THIS DATA INTO DISK FILE 12.
C
      TDSKWRK = DSKWRK
      DSKWRK = .FALSE.
C
      CALL SEQOPN(NFT12,'CIVECTR','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQREW(NFT12)
C
      DO IST=1,NSTATE
         LPOINT = LVEC + (IST-1)*NCFG
CSPW     IF(MASWRK) CALL STFASE(X(LPOINT),NCFG,NCFG,1)
         CALL SQWRIT(NFT12,X(LPOINT),NCFG)
      END DO
      CALL SEQREW(NFT12)
      DSKWRK = TDSKWRK
C
      CALL RETFM(NEED)
      RETURN
C
 9000 FORMAT(1X,'READING',I5,' STATES EXPRESSED IN TERMS OF',I10,1X,A4)
 9010 FORMAT(I10,1X,A4,1X,I4)
 9030 FORMAT(5E14.6)
 9040 FORMAT(1X,'ERROR, THERE IS A MISMATCH IN THE $CISVEC GROUP:'/
     *       5X,'EXPECTING ',I10,1X,A4,' FOUND',I10,1X,A4/
     *       5X,'EXPECTING ',I4,' STATES, FOUND',I4,' STATES')
      END
C*MODULE CISGRD  *DECK TRDENM
      SUBROUTINE TRDENM(C,X,A,T,TOCC,TVIR,DINTX,DINTY,DINTZ,
     &                 EIGEV,MULT,NCFG,NBF,NBF2,NOC,NVIR,NCOR,NOCA,
     &                 NSTATE,ESCF,MASWRK,IW,DMX,DMY,DMZ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL MASWRK
C
      DIMENSION C(NBF,NBF),X(NBF,NBF),A(NCFG,*),T(NBF,NBF)
      DIMENSION TOCC(NOC,NOC),TVIR(NVIR,NVIR)
      DIMENSION DINTX(NBF2),DINTY(NBF2),DINTZ(NBF2)
      DIMENSION EIGEV(NSTATE),MULT(NSTATE)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00,
     *           THREE=3.0D+00, EIGHT=8.0D+00, PT5=0.5D+00)
C
C     CONVERSION FACTORS FROM NATIONAL INSTITUTE OF STANDARDS AND TECH.
C
      PARAMETER (PLANCK=6.6260693D-27, CLIGHT=2.99792458D+10,
     *           BOHR=5.291772108D-09, ECHARG=4.803204401D-10,
     *           EV=27.2113845D+00, DEBYE=2.541766D+00)
C
C     --- CALCULATE CIS TRANSITION DIPOLE MOMENTS ---
C
      PI = ACOS(-ONE)
      ERG = ECHARG*ECHARG/BOHR
C
      IF(MASWRK) WRITE(IW,8000)
      IF(MASWRK) WRITE(IW,8100) DMX,DMY,DMZ
C
      DO 50 II=0,NSTATE
C
         IF(II.EQ.0) THEN
            FACTOR = SQRT(TWO)
         ELSE
            FACTOR = ONE
         END IF
C
         DO 60 JJ=II,NSTATE
C
            IF(JJ.EQ.0) GOTO 60
C
            IF(II.EQ.0) THEN
C
C           --- IF ONE OF THE STATES GROUND STATE SIMPLY TRANSFORM
C               CIS COEFICENTS OF OTHER STATE TO AO BASIS TO FORM T ---
C               (NOTE IN THIS CASE T MULTIPLIED BY ROOT TWO BELOW)
C
               CALL MRARTR(A(1,JJ),NVIR,NVIR,NOCA,C(1,NCOR+1),
     &                     NBF,NBF,X,NBF)
               CALL MRARBR(C(1,NOC+1),NBF,NBF,NVIR,X,NBF,NBF,T,NBF)
            ELSE
C
C       --- OTHERWISE FORM CROSS DENSITY MATRIX AND THEN
C                  TRANFORM TO AO BASIS TO FORM T ---
C
               CALL DGEMM('T','N',NOCA,NOCA,NVIR,-ONE,A(1,II),NVIR,
     &                    A(1,JJ),NVIR,ZERO,TOCC,NOCA)
               CALL MRARTR(A(1,II),NVIR,NVIR,NOCA,A(1,JJ),NVIR,
     &                     NVIR,TVIR,NVIR)
C
               CALL CISAOT(TOCC,TVIR,C,X,T,NBF,NCOR,NOCA,NOC,NVIR)
            END IF
C
            TDMX = ZERO
            TDMY = ZERO
            TDMZ = ZERO
C
            DO 100 I=1,NBF
               DO 110 J=1,I
                  TDEN = (T(I,J) + T(J,I))*FACTOR
                  IF(I.EQ.J) TDEN = TDEN * PT5
                  IJ = (I*I-I)/2 + J
                  TDMX = TDMX + TDEN*DINTX(IJ)
                  TDMY = TDMY + TDEN*DINTY(IJ)
                  TDMZ = TDMZ + TDEN*DINTZ(IJ)
  110          CONTINUE
  100       CONTINUE
C
            IF(II.EQ.JJ) THEN
               TDMX = -TDMX + DMX/DEBYE
               TDMY = -TDMY + DMY/DEBYE
               TDMZ = -TDMZ + DMZ/DEBYE
            END IF
C
C           --- PREPARE FINAL LENGTH FORM RESULTS ---
C
            IF(II.EQ.0) THEN
               MULST1 = 1
               ETOT1 =  ESCF
               DELENG = EIGEV(JJ)
            ELSE
               MULST1 = 2 * MULT(II) + 1
               ETOT1 = EIGEV(II) + ESCF
               DELENG = EIGEV(JJ) - EIGEV(II)
            END IF
C
            MULST2 = 2 * MULT(JJ) + 1
            ETOT2 = EIGEV(JJ) + ESCF
C
            EXC = EV*DELENG
            FREQ = DELENG*ERG/PLANCK
            OMEGA = FREQ/CLIGHT
C
            DIP = SQRT(TDMX*TDMX + TDMY*TDMY + TDMZ*TDMZ)
            WMX = TDMX*DEBYE
            WMY = TDMY*DEBYE
            WMZ = TDMZ*DEBYE
            WIP = DIP*DEBYE
C
C           --- OSCILLATOR STRENGTH AND EINSTEIN COEFFS ---
C
C           THE FORMULA FOR F ASSUMES DEL-E,<R> IN ATOMIC UNITS.
C
C                 8 * PI**3
C           B = -------------  * <R>**2
C                3 * H**2 * C
C
C                8 * PI * H * FREQ**3
C           A = ---------------------- * B
C                      C**2
C
C           THE ABOVE FORMULAE ASSUME G=1 FOR BOTH STATES.
C
            F = TWO * DELENG * DIP * DIP/THREE
            FACT = (DEBYE * 1.0D-18 / PLANCK)
            FACT = EIGHT * PI * PI * PI * FACT * FACT
            FACT = FACT/(THREE*CLIGHT)
            B = FACT * DIP * DIP
            AA = EIGHT * PI * DELENG * ERG * OMEGA * OMEGA * B
C
            IF(MASWRK) THEN
               IF(II.EQ.0) THEN
                  IF(MULST1.EQ.MULST2) THEN
                     WRITE(IW,9000) JJ
                     WRITE(IW,9010) MULST1,MULST2,ETOT1,ETOT2,
     &                              FREQ,OMEGA,EXC
                     WRITE(IW,9200) TDMX,TDMY,TDMZ,DIP,WMX,WMY,WMZ,WIP
                  ELSE
                     WRITE(IW,9001) JJ,ETOT2,MULST2,FREQ,OMEGA,EXC
                  ENDIF
               ELSE IF(II.EQ.JJ) THEN
                  WRITE(IW,9015) II
                  WRITE(IW,9016) MULST1,ETOT1
                  WRITE(IW,9201) TDMX,TDMY,TDMZ,DIP,WMX,WMY,WMZ,WIP
               ELSE
                  IF(MULST1.EQ.MULST2) THEN
                     WRITE(IW,9020) II,JJ
                     WRITE(IW,9030) MULST1,MULST2,ETOT1,ETOT2,
     &                              FREQ,OMEGA,EXC
                     WRITE(IW,9200) TDMX,TDMY,TDMZ,DIP,WMX,WMY,WMZ,WIP
                  ELSE
                     WRITE(IW,9021) II,MULST1,JJ,MULST2
                  ENDIF
               END IF
               IF(II.NE.JJ  .AND.  MULST1.EQ.MULST2)
     *               WRITE(IW,9210) F,AA,B
            END IF
C
   60    CONTINUE
   50 CONTINUE
      RETURN
C
 8000 FORMAT(/2X,69(1H-),/,20X,'CIS TRANSITION DIPOLE MOMENTS AND',
     &                   /,20X,'EXPECTATION VALUES OF DIPOLE MOMENTS',
     &       /2X,69(1H-))
 8100 FORMAT(/1X,'GROUND STATE (SCF) DIPOLE=',3F13.6,' DEBYE')
 9000 FORMAT(/1X,'TRANSITION FROM THE GROUND STATE TO EXCITED',
     &       ' STATE' ,I3)
 9001 FORMAT(/1X,'TRANSITION FROM THE GROUND STATE TO EXCITED STATE',
     *            I3,' IS OPTICALLY INACTIVE,'/
     *        1X,'STATE ENERGY =  ',F20.10,' MULT=',I3/
     *        1X,'EXCITATION ENERGY =',1P,E12.4,0P,' [1/SEC] =',F12.2,
     *          ' [1/CM] =',F10.2,' [EV]')
 9010 FORMAT(/1X,'STATE MULTIPLICITIES = ',2I3,
     *      /1X,'STATE ENERGIES =  ',2F20.10,
     *      /1X,'EXCITATION ENERGY =',1P,E12.4,0P,' [1/SEC] =',F12.2,
     *          ' [1/CM] =',F10.2,' [EV]')
 9015 FORMAT(/1X,'EXPECTATION VALUE DIPOLE MOMENT FOR EXCITED',
     &       ' STATE' ,I3)
 9016 FORMAT(/1X,'STATE MULTIPLICITY = ',I3,
     *      /1X,'STATE ENERGY =  ',F20.10)
 9020 FORMAT(/1X,'TRANSITION BETWEEN EXCITED STATES ',I3,' AND',I3)
 9021 FORMAT(/1X,'TRANSITION FROM STATE',I3,', MULT=',I2,
     *           ' TO STATE',I3,', MULT=',I2,
     *           ' IS OPTICALLY INACTIVE.')
 9030 FORMAT(/1X,'STATE MULTIPLICITIES=',2I3,
     *      /1X,'STATE ENERGIES =  ',2F20.10,
     *      /1X,'TRANSITION ENERGY =',1P,E12.4,0P,' [1/SEC] =',F12.2,
     *          ' [1/CM] =',F12.2,' [EV]')
 9200 FORMAT(1X,27X,'X ',10X,'Y ',10X,'Z ',10X,
     *              'NORM'/
     *       1X,'TRANSITION DIPOLE =',4F12.6,' E*BOHR'/
     *       1X,'TRANSITION DIPOLE =',4F12.6,' DEBYE')
 9201 FORMAT(1X,27X,'X ',10X,'Y ',10X,'Z ',10X,
     *              'NORM'/
     *       1X,'     STATE DIPOLE =',4F12.6,' E*BOHR'/
     *       1X,'     STATE DIPOLE =',4F12.6,' DEBYE')
 9210 FORMAT(1X,'OSCILLATOR STRENGTH =',F12.6/
     *       1X,'EINSTEIN COEFFICIENTS: A=',1P,E12.4,
     *          ' 1/SEC; B=',E12.4,' SEC/G')
      END
C
C*MODULE CISGRD  *DECK CISAOT
      SUBROUTINE CISAOT(PIJ,PAB,C,X,T,NBF,NCOR,NOCA,NOC,NVIR)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ONE=1.0D+00)
C
      DIMENSION PIJ(NOC,NOC),PAB(NVIR,NVIR)
      DIMENSION C(NBF,*),X(NBF,*),T(NBF,*)
C
C     --- TRANSFORM OCC-OCC AND VIR-VIR BLOCKS TO AO BASIS ---
C
      CALL MRARTR(PIJ,NOCA,NOCA,NOCA,C(1,NCOR+1),NBF,NBF,X,NBF)
      CALL MRARBR(C(1,NCOR+1),NBF,NBF,NOCA,X,NBF,NBF,T,NBF)
C
      CALL MRARTR(PAB,NVIR,NVIR,NVIR,C(1,NOC+1),NBF,NBF,X,NBF)
      CALL DGEMM('N','N',NBF,NBF,NVIR,ONE,C(1,NOC+1),NBF,X,NBF,ONE,
     &           T,NBF)
C
      RETURN
      END
C
