C 21 APR 08 - SHY - CHANGES FOR PARALLEL RUNS, LEBEDEV GRIDS
C 11 APR 08 - MWS - UTDOSCALC: DO DIPOLE INTEGRALS AT ALL GEOMETRIES
C  4 MAR 08 - SHY - NEW MODULE FOR SPIN UNRESTRICTED TD-DFT
C
C*MODULE UTDDFT  *DECK UTDDFTCALC
      SUBROUTINE UTDDFTCALC
C     ******************************************************************
C     PAPERS:
C     (1) S. HIRATA & M. HEAD-GORDON, CHEM. PHYS. LETT. 302, 375-382
C         (1999)
C     (2) R. BAUERNSCHMITT & R. AHLRICHS, CHEM. PHYS. LETT. 256, 454-464
C         (1996)
C     (3) R. BAUERNSCHMITT, M. HASER, O. TREUTLER, & R. AHLRICHS,
C         CHEM. PHYS. LETT. 264, 573-578 (1997)
C     (4) R. E. STRATMANN AND G.E. SCUSERIA, J. CHEM. PHYS. 108,
C        8218-8224 (1998).
C     COMMENT:
C     PAPER (1) : TDDFT FOR RADICALS (SPIN-UNRESTRICTED TDDFT)
C     PAPER (2) : TDDFT FOR CLOSE SHELLS, SEE APPENDIX FOR SINGLET AND
C                 TRIPLET EXCITATION
C     PAPER (3) : SEE 2. IMPLEMENTATION SINCE IT HAS DISCUSSED THE
C                 EFFICIENT PROCEDURE FOR TDDFT.
C     ******************************************************************
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXRT=100)
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00,TWO=2.0D+00)
C
      LOGICAL TRIPLET,SG1,SG1T,GOPARR,DSKWRK,MASWRK,PACK2E
C
      CHARACTER*20 METHOD
C
      PARAMETER (MXGRID=10)
C
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INFGRD/ RHOMIN,ILENG,MAXGRD
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /INFOTD/ TDM(3,MXRT),CNVTOL,TRIPLET,MOCC(MXRT),MVIR(MXRT),
     *                NRADT,NTHET,NPHIT,NSTAT,NTRIAL,MAXVEC,NTHST,
     *                IRECTD,NLEBT,SG1T
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA CHECK/8HCHECK   /
      DATA ROHF /8HROHF    /
      DATA UHF  /8HUHF     /
C
      DATA MXITER/100/
C
C     ----- DRIVER FOR OPEN SHELL + TDDFT EXCITATION ENERGIES -----
C
C        OF COURSE, MANY OF THE EXCITED STATES ARE NOT TOTALLY SYMMETRIC
C
      CALL SYMOFF
C
      CALL VCLR(TDM,1,3*MXRT)
C
C     --- DEFINE NOCC & NVIR
C
      LX   = NQMT
      NOCA = NA
      NOCB = NB
      L0   = NINTMX
      L1   = NUM
      L2   = L1*(L1+1)/2
      L3   = L1*L1
      L7A  = NOCA*(LX-NOCA)
      L7B  = NOCB*(LX-NOCB)
      L7C  = L7A+L7B
      L7   = MAX(L7A,L7B)
      CSCAL = ONE
      IPF  = 1
      IMF  = 0
      IF (DFTTYP(3).NE.ZERO) IMF=1
C     --- DEFINE MXVEC & NDSR & NVEC:NUMBER OF INITIAL VECORS
      MXVEC=MAXVEC
      NDSR=NSTAT
      NVEC=NTRIAL
C     --CHECK
      IF(NDSR .GT.L7)   NDSR=L7
      MXVEC=MAXVEC*NDSR
      IF(MXVEC.GT.L7)   MXVEC=L7
C     --LAST CHECK
      IF(NVEC .GT.MXVEC)NVEC=MXVEC
      IF(NVEC .LT.NDSR) NVEC=NDSR
C
      NMAX=MAX(2*NDSR,NTRIAL)
      MXVEC2=MXVEC*(2*MXVEC+1)  !(2*MXVEC*(2*MXVEC+1)/2)
      MXVEC3=MXVEC*MXVEC
      IST=1
      IEND=NVEC
C     == MEMORY
      CALL GOTFM(NGOTMX)
      CALL VALFM(LOADFM)
C
      LVA   = LOADFM+1
      LVB   = LVA   + L1*LX     ! VA(L1,LX): MO (A) : K70 IN TTDFT.SRC
      LEA   = LVB   + L1*LX     ! VB(L1,LX): MO (B)
      LEB   = LEA   + LX        ! EA(LX)   : EIGENVALUES FOR (A) ORBS
      LINXOV= LEB   + LX        ! EB(LX)   : EIGENVALUES FOR (B) ORBS
      LXM   = LINXOV+ L7C*2     ! INXOVA(2,L7C)
      LIM   = LXM   + L7C       ! XM(L7C)
      LBVA  = LIM   + L7C       ! IM(L7C)
      LBVB  = LBVA  + L7*MXVEC  ! BVA(L7,MXVEC)
      LWRK1 = LBVB  + L7*MXVEC  ! BVB(L7,MXVEC)
      LWRK2 = LWRK1 + L3        ! WRK1(L1,L1)
      LWRK3 = LWRK2 + L3        ! WRK2(L1,L1)
      L2EA  = LWRK3 + L3        ! WRK3(L1,L1)
      L2EB  = L2EA  + L3*NMAX   ! L2EA(L1,L1,NV)
      LBUF  = L2EB  + L3*NMAX   ! L2EB(L1,L1,NV)
      LIBUF = LBUF  + L0        ! XX(L0) IN TD2E
      LAPB  = LIBUF + L0*LABSIZ ! IX(L0*LABSIZ) IN TD2E
      LAMB  = LAPB  + L2*NMAX*2 ! APB(L2,2,NV)   := A+B FOR (A,A),(B,B)
      LPMOA = LAMB  + L3*NMAX*2 ! AMB(L1,L1,2,NV):= A-B FOR (A,A),(B,B)
      LPMOB = LPMOA + L7*MXVEC  ! PMOA(L7,NV) := A+B FOR ALPHA
      LMMOA = LPMOB + L7*MXVEC  !
      LMMOB = LMMOA + L7*MXVEC  ! MMOA(L7,NV) := A-B FOR ALPHA
      LRMAPB= LMMOB + L7*MXVEC  ! MMOB(L7,NV) := A-B FOR BETA
      LRMAMB= LRMAPB+ MXVEC3    ! XM_P(NV,NV):= REDUCED [A+B]
      LRMEE = LRMAMB+ MXVEC3    ! XM_M(NV,NV):= REDUCED [A-B]
      LRMVR = LRMEE + MXVEC     ! RMEE(MXVEC): EIGENVAL FOR REDUCED MTRX
      LRMVL = LRMVR + MXVEC3    ! RMVR(MXVEC,MXVEC): VR FOR REDUCED MTRX
      LRMSR2= LRMVL + MXVEC3    ! RMVL(MXVEC,MXVEC): VL FOR REDUCED MTRX
      LRMSRC= LRMSR2+ MXVEC*8   ! RMSC2(MXVEC*8) : VL FOR REDUCED MATRIX
      LIWRK = LRMSRC+ MXVEC2    ! RMSCR(MXVEC*8) : VL FOR REDUCED MATRIX
      LIIWRK= LIWRK + L7        ! RIWRK(L7) :
      LVRO  = LIIWRK+ MXVEC     ! RIIWRK(MXVEC) :
      LVLO  = LVRO  + L7*NDSR*2 ! VRO(L7,NDSR,2) :
      LABP  = LVLO  + L7*NDSR*2 ! VLO(L7,NDSR,2) :
      LABM  = LABP  + L7*NDSR*2 ! ABP_R(L7,NDSR,2) :
      LQA   = LABM  + L7*NDSR*2 ! ABM_L(L7,NDSR,2) :
      LQB   = LQA   + L7*NDSR*2 ! QA (L7,NDSR,2) :
      LABSYM= LQB   + L7*NDSR*2 ! QB (L7,NDSR,2) :
      LTXYZ = LABSYM+ NDSR      ! LABSYM(NDSR)
      LOS   = LTXYZ + 3*NDSR    ! TXYZ(3,NDSR)
      LAX   = LOS   + NDSR      ! OS (NDSR)
      LAY   = LAX   + L2
      LAZ   = LAY   + L2
      LIRMOA= LAZ   + L2
      LIRMOB= LIRMOA+ LX
      LAST  = LIRMOB+ LX
C     MEMORY FOR DFT GRID
      IF (NDFTFG.EQ.1) THEN
         LGRD  = LAST
         LWGT  = LGRD + MAXGRD*3 ! XYZGRD(MAXGRD*3)
         LDCH  = LWGT + MAXGRD   ! XYZWGT(MAXGRD)
         LRHOI = LDCH + 4*L1*ILENG ! XYZDCH
         LTRAI = LRHOI+ 8*MAXGRD ! DRHOI_(MAXGRD,4,2)
         LCOEF = LTRAI+ 8*ILENG  ! TRAI (4,2,ILENG)
         LEX   = LCOEF+ 16*ILENG ! COEF (4,4,ILENG)
         LEX0  = LEX  + 10*ILENG ! EX (ILENG,10), IN CASE OF 2ND DERIVS
         LEC   = LEX0 + ILENG    ! EX0(ILENG)
         LEC0  = LEC  + 20*ILENG ! EC (ILENG,20), IN CASE OF 2ND DERIVS
         LAOMAX= LEC0 + ILENG    ! EC0(ILENG)
         LVHA  = LAOMAX + L1    ! AOMAX(L1)
         LVHB  = LVHA + LX*L1   ! VHA(LX,L1)
         LGMO  = LVHB + LX*L1   ! VHB(LX,L1)
         LVALGA= LGMO + L1*8    ! GMO (L1,4,2) FOR ALPHA AND BETA
         LFXC  = LVALGA+ L1*2   ! VALGA (L1,2) FOR ALPHA AND BETA
         LAST  = LFXC + L2*2    ! FXC(L2,2)
      END IF
C
      NEED = LAST-LOADFM-1
      IF(MASWRK) THEN
         IF (SCFTYP.EQ.ROHF) THEN
            WRITE(UNIT=METHOD,FMT='(''TD-DFT/ROHF-'',A8)') DFTYPE
         ELSE IF (SCFTYP.EQ.UHF) THEN
            WRITE(UNIT=METHOD,FMT='(''TD-DFT/UHF -'',A8)') DFTYPE
         ELSE
            WRITE(UNIT=METHOD,FMT='(''ERROR IN TD-DFT'')')
            GO TO 1000
         END IF
         WRITE(IW,9010) METHOD,NEED
      END IF
      CALL FLSHBF(IW)
      CALL GETFM(NEED)
C
      IF(EXETYP.EQ.CHECK) GO TO 800
C     --- READ MO (LVA)
      CALL DAREAD(IDAF,IODA,X(LVA),L1*LX,15,0)
      IF(SCFTYP.EQ.ROHF) THEN
         CALL DCOPY(L1*LX,X(LVA),1,X(LVB),1) ! FOR ROHF, LVB=LVA
      ELSE
         CALL DAREAD(IDAF,IODA,X(LVB),L1*LX,19,0) !
      END IF
C
C     --- GRID PART
C
      IF (NDFTFG.EQ.1) THEN
C     --- CREATE DRHOI_ & XC DERIVATIVES
         NPTGRD=MAXGRD
         CALL UTDDFTSET(X(LGRD),X(LWGT),X(LDCH),X(LVA),X(LVB),
     *        X(LRHOI),X(LAOMAX),X(LGMO),ILENG,NPTGRD,L1)
      END IF
C
C     --- INITIAL TRIAL VECTORS BVEC(L7,NVEC)
C     READ EIGENVALUES FOR ALPHA ORBITAL
      CALL DAREAD(IDAF,IODA,X(LEA),LX,17,0)
      IF(SCFTYP.EQ.ROHF) THEN
         CALL DCOPY(LX,X(LEA),1,X(LEB),1) ! FOR ROHF, LEB=LEA
      ELSE
         CALL DAREAD(IDAF,IODA,X(LEB),LX,21,0) ! FOR UHF, LEA!=LEB
      END IF
C
      CALL UTDINIVEC(X(LEA),X(LEB),X(LINXOV),X(LXM),X(LIM),
     *               X(LBVA),X(LBVB),LX,L7,L7C,NOCA,NOCB,NVEC,
     *               MXVEC,X(LIIWRK))
C
      IF(MASWRK) WRITE(IW,*)'INITIAL VECTORS / MXVEC ',NVEC,'/',MXVEC
C     --NEW IEND
      IEND=NVEC
C
C     ---- TDDFT LOOP START ----
C
      ITER=0
 10   CONTINUE
      ITER=ITER+1
      IF(MASWRK) WRITE(IW,*) 'ITERATION # ', ITER
      IF(ITER.GT.MXITER) THEN
         IF(MASWRK) THEN
            WRITE(IW,*) '  .. OVER MAXIMUM ITERATION NUMBER .. '
            WRITE(IW,*) ' ITER, MXITER =', ITER, MXITER
         END IF
         ICK = 1
         GOTO 20
      END IF
C ********************************************************************
C
C     --- MATRIX-VECTOR PRODUCTS
C         (A+B)_B & (A-B)_B
C
      NV=IEND-IST+1
      IF(MASWRK) WRITE(IW,9020) IST,IEND,NV
      CALL FLSHBF(IW)
C
C     -- AO BASIS
C
      IV = 0
      DO IVEC=IST, IEND
         IV=IV+1
C     TRANSFORM MO TO AO
C     FOR ALPHA:
         CALL IATOGENTD(X(LBVA),X(LWRK1),LX,L7,NOCA,IVEC,MXVEC)
C    (1) LVA*BVA*LVA^T
         CALL DGEMM('N','N',L1,LX,LX,ONE,X(LVA),L1,X(LWRK1),LX,ZERO,
     *        X(LWRK2),L1)
         CALL DGEMM('N','T',L1,L1,LX,ONE,X(LWRK2),L1,X(LVA),L1,ZERO,
     *        X(L2EA+(IV-1)*L1*L1),L1)
C
C     FOR BETA
C
         CALL IATOGENTD(X(LBVB),X(LWRK1),LX,L7,NOCB,IVEC,MXVEC)
C     (1) LVB*BVB*LVB^T
         CALL DGEMM('N','N',L1,LX,LX,ONE,X(LVB),L1,X(LWRK1),LX,ZERO,
     *        X(LWRK2),L1)
         CALL DGEMM('N','T',L1,L1,LX,ONE,X(LWRK2),L1,X(LVB),L1,ZERO,
     *        X(L2EB+(IV-1)*L1*L1),L1)
      END DO
C
C     - 2E : A+B AND A-B
C
C     [A+B](A) AND [A-B](A)
C     [A+B](A,A) => 2 (IA(A)|MU VU) - C_HF (I MU | A VU)
C                                   - C_HF (I VU | A MU)
C     [A+B](A,B) => 2 (IA(A)|MU VU)
C     [A+B](B,A) => 2 (IA(B)|MU VU)
C     [A+B](B,B) => 2 (IA(B)|MU VU) - C_HF (I MU | A VU)
C                                   - C_HF (I VU | A MU)
C     [A-B](A,A) => - C_HF (I MU | A VU) + C_HF (I VU| A MU)
C     [A-B](B,B) => - C_HF (I MU | A VU) + C_HF (I VU| A MU)
      CALL VCLR(X(LAPB),1,2*L2*NMAX)
      CALL VCLR(X(LAMB),1,2*L3*NMAX)
      LAPBA=LAPB
      LAPBB=LAPB+L2*NMAX
      LAMBA=LAMB
      LAMBB=LAMB+L3*NMAX
      CALL UTD2E(X(L2EA),X(L2EB),X(LAPBA),X(LAPBB),
     *     X(LAMBA),X(LAMBB),X(LBUF),X(LIBUF),L1,CSCAL,IMF,IPF,NV)
C
C     --- DFT CONTRIBUTION
C
      IV = 0
      DO IVEC=IST,IEND
         IV=IV+1
         LAPB1 = LAPBA+(IV-1)*L2
         LAPB2 = LAPBB+(IV-1)*L2
         IF (NDFTFG.EQ.1) THEN
C     -- DFT FUNCTIONAL PART
C     -- HALF-TRANSFORMING
C     : IN ROUTINE UTDFXCP, LVA (OR LVB) WILL BE MULTIPLIED TO
C     BVA*LVA^T (OR BVB*LVB^T)
C     COMMENT:========================================
C     (2) TRANSFORM Z_IA (MO) INTO Z_IA(AO).
C     ================================================
            CALL IATOGENTD(X(LBVA),X(LWRK1),LX,L7,NOCA,IVEC,MXVEC)
C     (1) BVA*LVA^T
            CALL DGEMM('N', 'T', LX, L1, LX, ONE, X(LWRK1), LX,
     *           X(LVA), L1, ZERO, X(LVHA), LX)
C
            CALL IATOGENTD(X(LBVB),X(LWRK1),LX,L7,NOCB,IVEC,MXVEC)
C     (1) BVB*LVB^T
            CALL DGEMM('N', 'T', LX, L1, LX, ONE, X(LWRK1), LX,
     *           X(LVB), L1, ZERO, X(LVHB), LX)
C     (3) PERFORM THE MULTIPLICATION IN THE AB BASIS
C         FXC(MU,VU,1) = (MU VU(A)|W(A,A)|JB (A)) +
C                        (MU VU(A)|W(A,B)|JB (B))
C         FXC(MU,VU,2) = (MU VU(B)|W(B,A)|JB (A)) +
C                        (MU VU(B)|W(B,B)|JB (B))
            CALL VCLR(X(LFXC),1,L2*2)
            CALL UTDFXCP(X(LFXC),RHO,X(LGRD),X(LWGT),X(LDCH),
     *           X(LVHA),X(LVHB),X(LVA),X(LVB),X(LRHOI),X(LTRAI),
     *           X(LCOEF),X(LEX),X(LEC),X(LEX0),X(LEC0),
     *           X(LAOMAX),X(LGMO),X(LVALGA),
     *           ILENG,NPTGRD,L1,L2,LX,2)
C     NOW WE HAVE FXC(L2,4) !1=>AA,2=>AB,3=>BA,4=>BB
C     NOTE, [A+B]_XC = (MU VU|W|JB) + (MU VU|W|BJ)
C                    = 2 (MU VU|W|JB) IN CASE OF REAL WAVEFUNCTIONS.
            LFXC1 = LFXC
            LFXC2 = LFXC+L2
            CALL DAXPY(L2,TWO,X(LFXC1),1,X(LAPB1),1)
            CALL DAXPY(L2,TWO,X(LFXC2),1,X(LAPB2),1)
         END IF
C     =====================================
C     (4) TRANSFORM BACK INTO THE MO BASIS:
C     =====================================
C     [A+B](IA(A)) = [A+B](IA(A),JB(A)) + [A+B](IA(A),JB(B)) := PMOA
C     [A+B](IA(B)) = [A+B](IA(B),JB(A)) + [A+B](IA(B),JB(B)) := PMOB
         CALL UMNTOIAP(X(LAPB1),X(LAPB2),X(LPMOA),X(LPMOB),
     *                 X(LVA),X(LVB),X(LWRK1),X(LWRK2),
     *                 L1,L2,LX,NOCA,NOCB,L7,IVEC)
         CALL ESUM(X(LEA),X(LPMOA),X(LBVA),LX,L7,NOCA,IVEC)
         CALL ESUM(X(LEB),X(LPMOB),X(LBVB),LX,L7,NOCB,IVEC)
C     PRODUCT A-B : FROM UTD2E
C     [A-B](IA(A)) = [A-B](IA(A),JB(A)) := MMOA
C     [A-B](IA(B)) = [A-B](IA(B),JB(B)) := MMOB
         LAMB1 = LAMBA+(IV-1)*L3
         LAMB2 = LAMBB+(IV-1)*L3
         CALL UMNTOIAM(X(LAMB1),X(LAMB2),X(LMMOA),X(LMMOB),
     *                 X(LVA),X(LVB),X(LWRK1),L1,LX,NOCA,NOCB,L7,IVEC)
         CALL ESUM(X(LEA),X(LMMOA),X(LBVA),LX,L7,NOCA,IVEC)
         CALL ESUM(X(LEB),X(LMMOB),X(LBVB),LX,L7,NOCB,IVEC)
      END DO
C     *************************************
C     --- SMALL REDUCED MATRICES A+B, A-B
C     *************************************
C     REDUCED MATRICES
C     [A+B](P,Q) = Z_IA(A)^P * [A+B]_IA(A)^Q + Z_IA(B)^P * [A+B]_IA(B)^Q
C     [A-B](P,Q) = Z_IA(A)^P * [A-B]_IA(A)^Q + Z_IA(B)^P * [A-B]_IA(B)^Q
C     FOR [A+B](P,Q)(A) = Z_IA(A)^P * [A+B]_IA(A)^Q AND
C         [A-B](P,Q)(A) = Z_IA(A)^P * [A-B]_IA(A)^Q
C     FOR [A+B](P,Q)(B) = Z_IA(B)^P * [A+B]_IA(B)^Q AND
C         [A-B](P,Q)(B) = Z_IA(B)^P * [A-B]_IA(B)^Q
C
      CALL URPAREDMS(X(LBVA),X(LBVB),X(LPMOA),X(LMMOA),
     *               X(LPMOB),X(LMMOB),X(LRMAPB),X(LRMAMB),
     *               NVEC,NOCA,NOCB,LX,L7)
C
C     --- FORM REDUCED RPA MATRIX AND DIAGONALIZE IT
C     FOR ALPHA ORBITAL
C
      CALL RPAEIG(X(LRMEE), X(LRMVL), X(LRMVR), X(LRMAPB), X(LRMAMB),
     *            X(LRMSR2),X(LRMSRC),X(LIWRK), NVEC)
C
C     --- IMPOSING BIORTHOGONALITY CONDITION ON V_L & V_R
C
      CALL RPAVNORM (X(LRMVR), X(LRMVL), NVEC) ! FOR ALPHA+BETA
C
C     --- PRINT CURRENT EXCITATION ENERGIES
C
      CALL URPAECHK(X(LRMEE),X(LIWRK),X(LIIWRK),NVEC,NDSR,IMAX)
C
C     --- VRO(L7,*), VLO(L7,*), (A+B)VRO, (A-B)VLO
C         LVRO       LVLO       LABP      LABM
C     FOR ALPHA ORBITALS
C
      CALL URPAEXPNDV(X(LRMVR),X(LRMVL),X(LVRO),X(LVLO),
     *                X(LBVA),X(LBVB),NDSR,NVEC,L7A,L7B,L7)
      CALL URPAEXPNDV2(X(LRMVR),X(LRMVL),X(LABP),X(LABM),
     *                 X(LPMOA),X(LPMOB),X(LMMOA),X(LMMOB),
     *                 NDSR,NVEC,L7A,L7B,L7)
C
C     --- RESIDUAL VECTORS W AND CHECK CONVERGENCE
C         AND GET PERTURBED VECTORS Q IF REQUIRED
C
      CALL URPARESVEC(X(LQA),X(LQB),X(LABP),X(LABM),X(LVLO),X(LVRO),
     *                X(LRMEE),X(LXM),L7A,L7B,L7,NDSR,ICK,ERROR,
     *                CNVTOL,IMAX)
C
C     --- FOR CONVERGED CASE
C
      IF(ICK.EQ.0) GOTO 20
C
C     ORTHONORMALIZE Q(L7,NDSR*2) AND APPEND NEW B VECTORS
C
      CALL URPANEWB(L7,L7A,L7B,NDSR,MXVEC,X(LBVA),X(LBVB),X(LQA),X(LQB),
     *              NOVEC,NVEC,ICK,X(LIIWRK))
C
      IF(ICK.EQ.1) THEN
         IF(MASWRK) THEN
            WRITE(IW,*)'  ..SOMETHING IS WRONG.. NVEC.GT.MXVEC'
            WRITE(IW,*)'NVEC/MXVEC=',NVEC,'/',MXVEC
         ENDIF
         GOTO 20
      ENDIF
C
C     --- UPDATE IST AND IVEC
      IST=NOVEC+1
      IEND=NVEC
      GOTO 10
C
 20   CONTINUE
      IF (ICK.EQ.0) THEN
         IF(MASWRK) THEN
            WRITE(IW,*) ' '
            WRITE(IW,*) '   ------------------------------------------'
            WRITE(IW,9040) METHOD
            WRITE(IW,*) '   ------------------------------------------'
         ENDIF
      ELSE
         IF(MASWRK) THEN
            WRITE(IW,*)' '
            WRITE(IW,*)'  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
            WRITE(IW,9041) METHOD
            WRITE(IW,*)'  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
         ENDIF
      END IF
C
C     ---- TDDFT CALCULATION FINISHED ----
C     RECOVER POINT GROUP SYMMETRY SO THAT WE CAN LABEL THE STATES
C
      CALL SYMON
C
      IF(MASWRK) THEN
         WRITE(IW,9977) NVEC ,MXVEC
         WRITE(IW,9978) ITER, MXITER
         WRITE(IW,9979) L7
         WRITE(IW,9980) NUM
         IF(NDFTFG.EQ.1) WRITE(IW,2111) RHO
      ENDIF
C
C     --- OSCILLATOR STRENGTH
C
C     LWRK1->TDENA, LWRK2->TDENB,LWRK3->SCR
      CALL UTDOSCALC(X(LOS),X(LTXYZ),X(LWRK1),X(LWRK2),X(LWRK3),
     *     X(LVRO),X(LVLO),X(LRMEE),X(LVA),X(LVB),X(LAX),X(LAY),X(LAZ),
     *     L1,LX,L2,L7,L7A,L7C,NDSR,X(LINXOV))
C
C         DON'T PRINT RESULTS IF WE DIDN'T CONVERGE, TO EMPHASIZE THAT
C
      IF(ICK.EQ.1) THEN
         IF(MASWRK) WRITE(IW,*) 'FAILED IN TDDFT STEP. STOP'
      ELSE
C     READ MO SYMMETRY INFORMATION
         CALL DAREAD(IDAF,IODA,X(LIRMOA),LX,356,1)
         CALL DAREAD(IDAF,IODA,X(LIRMOB),LX,357,1)
C
         CALL UTDPRNT(X(LOS),X(LABSYM),X(LTXYZ),X(LINXOV),NDSR,
     *        X(LRMEE),X(LVRO),X(LIRMOA),X(LIRMOB),L7,L7A,L7C,LX,IMAX)
      ENDIF
      CALL FLSHBF(IW)
C
 800  CONTINUE
      CALL RETFM(NEED)
      IF(MASWRK) WRITE(IW,9030) METHOD
 1000 CONTINUE
      CALL TIMIT(1)
      RETURN
C
 2111 FORMAT(' NUMBER OF TOTAL ELECTRONS    =',F20.7)
 9010 FORMAT(1X,'MEMORY REQUIRED FOR ',A20,'ITERATIONS IS',I10,
     *          ' WORDS.')
 9020 FORMAT(1X,'   CURRENT VECTORS ',I4,' TO',I4,' (',I4,' VECTORS)')
 9030 FORMAT(1X,'..... DONE WITH ',A20,' EXCITATION ENERGIES .....')
 9040 FORMAT(1X,'   ', A20, ' CALCULATION CONVERGED ')
 9041 FORMAT(1X,'   ', A20, ' CALCULATION NOT CONVERGED ')
 9977 FORMAT(' NUMBER OF USED VECTORS / MAX =',I4,1X,'/',I4)
 9978 FORMAT(' NUMBER OF ITERATIONS   / MAX =',I4,1X,'/',I4)
 9979 FORMAT(' NUMBER OF SINGLE EXCITATIONS =',I6)
 9980 FORMAT(' NUMBER OF ATOMIC ORBITALS    =',I4)
      END
C
C*MODULE UTDDFT  *DECK UTDINIVEC
      SUBROUTINE UTDINIVEC(EA,EB,INXOV,XM,IM,BVECA,BVECB,LX,L7,L7C,
     *                     NOCA,NOCB,NVEC,MXVEC,IINDX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL DBUG,GOPARR,DSKWRK,MASWRK
      DIMENSION EA(LX),  EB(LX)
      DIMENSION XM(L7C)
      DIMENSION IM(L7C)
      DIMENSION BVECA(L7,*)
      DIMENSION BVECB(L7,*)
      DIMENSION INXOV(2,L7C)
      DIMENSION IINDX(*)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      PARAMETER (TOEV=27.21138386D+00)
      DATA MXINIVEC/50/
      DBUG=.FALSE.
C     -- MXINIVEC
      IF(MXINIVEC.GT.MXVEC) MXINIVEC=MXVEC
C     -- SET INXOV(2,L7C) & XM(L7C)
C     FOR ALPHA
      IJ=0
      DO J=NOCA+1,LX
         DO I=1,NOCA
            IJ=IJ+1
            INXOV(1,IJ)=I ! OCCUPIED ORBITAL
            INXOV(2,IJ)=J ! VIRTURAL ORBITAL
            XM(IJ)=EA(J)-EA(I)
            IM(IJ) =IJ    ! SAVE RPEVIOUS ORDERING
         ENDDO
      ENDDO
      L7A=IJ
C     FOR BETA ORBITALS
      DO J = NOCB+1,LX
         DO I=1,NOCB
            IJ=IJ+1
            INXOV(1,IJ)=I
            INXOV(2,IJ)=J
            XM(IJ)=EB(J)-EB(I)
            IM(IJ) =IJ    ! SAVE RPEVIOUS ORDERING
         END DO
      END DO
C     -- ORDERING XM(L7):XM(SMALL) =< XM(LARGE)
C     ** GET SMALLER DIAGONAL VALUES
      DO J=1,L7C-1
         DO I=J+1,L7C
            IF(XM(J).LE.XM(I)) GOTO 100
            XMJ=XM(J)
            IMJ=IM(J)
            XM(J)=XM(I)
            IM(J)=IM(I)
            XM(I)=XMJ
            IM(I)=IMJ
 100        CONTINUE
         END DO
      END DO
C     -- DBUG
      IF (DBUG) THEN
         WRITE(6,*) 'PRINT XM(L7) ORDERING'
         DO I=1,20
            WRITE(6,*) 'I,XM(I),IM(I)=',I,XM(I),IM(I)
         ENDDO
      ENDIF
C     -- SET XM(L7) AGAIN
      IJ=0
      DO J=NOCA+1,LX
         DO I=1,NOCA
            IJ=IJ+1
            XM(IJ)=EA(J)-EA(I)
         ENDDO
      ENDDO
C
      DO J=NOCB+1,LX
         DO I=1,NOCB
            IJ=IJ+1
            XM(IJ)=EB(J)-EB(I)
         ENDDO
      ENDDO
C
C     ----  GET INITIAL VECTORS
C     BVEC(L7,NVEC)
      CALL VCLR(BVECA,1,L7*MXVEC)
      CALL VCLR(BVECB,1,L7*MXVEC)
      DBUG=.TRUE.
      IF(DBUG.AND.MASWRK) THEN
         WRITE (6,'(14X,A6,1X,A5,6X,A5,2X,A10,4X,A4)')
     *        ' SPIN', 'NOCC', 'NVIR', ' (HARTREE) ', '(EV)'
      END IF
      DO K=1,NVEC
         IJ=IM(K)
         IF(IJ.LE.L7A) THEN
            IF(DBUG.AND.MASWRK)
     *           WRITE(6,'(A12,I2,A6,I6,A6,I5,F12.6,F10.4)')
     *           'TRIAL VEC.',K,'(A)',INXOV(1,IJ),'--->',
     *           INXOV(2,IJ),XM(IJ),XM(IJ)*TOEV
            BVECA(IJ,K) = 1.0D+00
            IINDX(K)=1
         ELSE
            IF(DBUG.AND.MASWRK)
     *           WRITE(6,'(A12,I2,A6,I6,A6,I5,F12.6,F10.4)')
     *           'TRIAL VEC.',K,'(B)',INXOV(1,IJ),'--->',
     *           INXOV(2,IJ),XM(IJ),XM(IJ)*TOEV
            BVECB(IJ-L7A,K) = 1.0D+00
            IINDX(K)=2
         END IF
      ENDDO
C
      END
C
C*MODULE UTDDFT  *DECK UMNTOIAP
      SUBROUTINE UMNTOIAP(PAOA,PAOB,PMOA,PMOB,VA,VB,QAO,SCR1,
     *                    L1,L2,LX,NOCA,NOCB,L7,IVEC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
C
      DIMENSION PAOA(L2),PAOB(L2)
      DIMENSION PMOA(L7,*)
      DIMENSION PMOB(L7,*)
      DIMENSION VA  (L1,LX)
      DIMENSION VB  (L1,LX)
      DIMENSION QAO (L1,L1)
      DIMENSION SCR1(*)!(L1,L1)
C
      CALL VCLR(PMOA(1,IVEC),1,L7)
      CALL VCLR(PMOB(1,IVEC),1,L7)
C
C     ==== FOR ALPHA ORBITAL, [A+B](IA(A)) =====
C     (1) PAO (L2) --> QAO(L1,L1)
      CALL CANTOGEN(PAOA,QAO,L1)
C
C     -- (M,N) TO (I,N)
C     SCRA = VA^T * QOA
C
      CALL DGEMM('T','N',NOCA,L1,L1,ONE,VA,L1,QAO,L1,ZERO,SCR1,NOCA)
C
      IJ=0
      DO J=NOCA+1,LX
         DO I=1,NOCA
            IJ=IJ+1
            DUMN=ZERO
            DO N=1,L1
               IN=(N-1)*NOCA+I
               DUMN=DUMN+SCR1(IN)*VA(N,J)
            ENDDO
            PMOA(IJ,IVEC)=DUMN
         ENDDO
      ENDDO
C
C     ==== FOR BETA ORBITAL, [A+B](IA(B)) =====
C     (1) PAO (L2) --> PAO(L1,L1)
C
      CALL CANTOGEN(PAOB,QAO,L1)
C
C     -- (M,N) TO (I,N)
C     SCRA = VA^T * QOA
C
      CALL DGEMM('T','N',NOCB,L1,L1,ONE,VB,L1,QAO,L1,ZERO,SCR1,NOCB)
C
      IJ=0
      DO J=NOCB+1,LX
         DO I=1,NOCB
            IJ=IJ+1
            DUMN=ZERO
            DO N=1,L1
               IN=(N-1)*NOCB+I
               DUMN=DUMN+SCR1(IN)*VB(N,J)
            ENDDO
            PMOB(IJ,IVEC)=DUMN
         ENDDO
      ENDDO
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK UMNTOIAM
      SUBROUTINE UMNTOIAM(PAOA,PAOB,MMOA,MMOB,VA,VB,SCR,
     *                    L1,LX,NOCA,NOCB,L7,IVEC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION MMOA, MMOB
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
      DIMENSION PAOA(L1,L1),PAOB(L1,L1)
      DIMENSION MMOA(L7,*),MMOB(L7,*)
      DIMENSION VA(L1,LX),VB(L1,LX)
      DIMENSION SCR(*)
C
C     ==== FOR ALPHA ORBITAL, [A-B](IA(A)) =====
C     (1) PAO (L1,L1,IV) --> PAO(L1,L1)
C
      CALL VCLR(MMOA(1,IVEC),1,L7)
      CALL VCLR(MMOB(1,IVEC),1,L7)
C
C     -- (M,N) TO (I,N)
C     SCR = VA^T * QAO
      CALL DGEMM('T','N',NOCA,L1,L1,ONE,VA,L1,PAOA,L1,ZERO
     *     ,SCR,NOCA)
C
C     -- (I,N) TO (IJ)
      IJ=0
      DO J=NOCA+1,LX
         DO I=1,NOCA
            IJ = IJ+1
            DUMN=ZERO
            DO N=1,L1
               IN=(N-1)*NOCA+I
               DUMN=DUMN+SCR(IN)*VA(N,J)
            ENDDO
            MMOA(IJ,IVEC)=DUMN
         ENDDO
      ENDDO
C
C
C     ==== FOR BETA ORBITAL, [A-B](IA(B)) =====
C     [A-B](IA(B)|MU,VU(B))
C     SCR = VB^T * QAO
      CALL DGEMM('T','N',NOCB,L1,L1,ONE,VB,L1,PAOB,L1,ZERO
     *     ,SCR,NOCB)
      IJ=0
      DO J=NOCB+1,LX
         DO I=1,NOCB
            IJ=IJ+1
            DUMN=ZERO
            DO N=1,L1
               IN=(N-1)*NOCB+I
               DUMN=DUMN+SCR(IN)*VB(N,J)
            ENDDO
            MMOB(IJ,IVEC)=DUMN
         ENDDO
      ENDDO
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK URPARESVEC
      SUBROUTINE URPARESVEC(QA,QB,W_L,W_R,V_L,V_R,EE,XM,
     *                      L7A,L7B,L7,NDSR,ICNVG,ERROR,TOL,IMAX)
C
C      CONSTRUCT RESIDUAL VECTORS AND CHECK CONVERGENCE
C     *INPUT W_L: ABP_R
C     *INPUT W_R: ABM_L
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      LOGICAL DBUG
C
      DIMENSION QA(L7,NDSR,2)
      DIMENSION QB(L7,NDSR,2)
      DIMENSION W_L(L7,NDSR,2)
      DIMENSION W_R(L7,NDSR,2)
      DIMENSION V_L(L7,NDSR,2)
      DIMENSION V_R(L7,NDSR,2)
      DIMENSION EE(*)
      DIMENSION XM(*)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      DBUG=.FALSE.
C
C     --- RESIDUAL VECTOR W_L & W_R
C
      DO IVEC=1,NDSR
C     --- ALPHA
         CALL DAXPY(L7A,-EE(IVEC),V_L(1,IVEC,1),1,W_L(1,IVEC,1),1)
         CALL DAXPY(L7A,-EE(IVEC),V_R(1,IVEC,1),1,W_R(1,IVEC,1),1)
C     --- BETA
         CALL DAXPY(L7B,-EE(IVEC),V_L(1,IVEC,2),1,W_L(1,IVEC,2),1)
         CALL DAXPY(L7B,-EE(IVEC),V_R(1,IVEC,2),1,W_R(1,IVEC,2),1)
      ENDDO
C
C     --NORMS OF W
C
      DO IVEC=1,NDSR
         CALL VCLR(QA(1,IVEC,1),1,L7)
         CALL VCLR(QA(1,IVEC,2),1,L7)
         CALL VCLR(QB(1,IVEC,1),1,L7)
         CALL VCLR(QB(1,IVEC,2),1,L7)
      END DO
C
      ERROR_MAX=0.0D+00
C
      DO IVEC=1,NDSR
         IF (IVEC.LE.IMAX) GOTO 100
C
C     --- GET NEW VECTORS Q
C     ALPHA
         ER_LA=DDOT(L7A,W_L(1,IVEC,1),1,W_L(1,IVEC,1),1)
         ER_RA=DDOT(L7A,W_R(1,IVEC,1),1,W_R(1,IVEC,1),1)
C
         ER_LB=DDOT(L7B,W_L(1,IVEC,2),1,W_L(1,IVEC,2),1)
         ER_RB=DDOT(L7B,W_R(1,IVEC,2),1,W_R(1,IVEC,2),1)
C
         ER_L =(ER_LA+ER_LB)
         ER_R =(ER_RA+ER_RB)
         ERR  =MAX(ER_L,ER_R)
C
         IF(DBUG) WRITE(IW,*) 'ERROR: IVEC/ERL/ERR=', IVEC, ER_L, ER_R
C
         IF(ERR.GT.1.0D+000) THEN
            WRITE(IW,*) 'LARGE ERROR DETECTED'
            WRITE(IW,*) 'STATE#=',IVEC
            WRITE(IW,*) 'ERROR LEFT/RIGHT=',ER_L,'/',ER_R
            ERR=0.0D+00
         ELSE IF(ERR.GT.ERROR_MAX)THEN
            ERROR_MAX=ERR
         END IF
C
         IF(ERR.GT.TOL) THEN
C     FOR ALPHA
            EEA=EE(IVEC)
            DO IJ=1,L7A
               DUM=1.0D+00/(EEA-XM(IJ)) !ABD(I))
               QA(IJ,IVEC,1)=DUM*W_R(IJ,IVEC,1)
               QA(IJ,IVEC,2)=DUM*W_L(IJ,IVEC,1)
            ENDDO
C     FOR BETA
            EEB=EE(IVEC)
            DO IJ=1,L7B
               DUM=1.0D+00/(EEB-XM(IJ+L7A)) !ABD(I))
               QB(IJ,IVEC,1)=DUM*W_R(IJ,IVEC,2)
               QB(IJ,IVEC,2)=DUM*W_L(IJ,IVEC,2)
            END DO
         ENDIF
 100     CONTINUE
      END DO
C
C     --CHECK CONVERGENCE OF TDDFT CALCULATION
C     ICNVG=0:CONVERGED
C
      IF(MASWRK) THEN
         WRITE(IW,1111) ERROR_MAX,TOL
         CALL FLSHBF(IW)
      ENDIF
C
      ERROR=ERROR_MAX
      IF (ERROR.LE.TOL) THEN
         ICNVG=0
      ELSE
         ICNVG=1
      ENDIF
C
      RETURN
 1111 FORMAT(1X,'   MAX ERROR =',1X,1P,E10.3,1X,'/',1P,E10.3)
      END
C
C*MODULE UTDDFT  *DECK URPANEWB
      SUBROUTINE URPANEWB(L7,L7A,L7B,NDSR,MXVEC,BVA,BVB,QA,QB,
     *     NOVEC,NVEC,ICK,IINDX)
C
C      ORTHONORMALIZE Q(NDIM,NDSR*2) AND APPEND TO BVEC
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL DBUG
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DIMENSION QA(L7,NDSR*2), QB(L7,NDSR*2)
      DIMENSION BVA(L7,*),     BVB(L7,*)
      DIMENSION IINDX(*)
      DATA THREB /1.0D-08/
      DBUG=.FALSE.
C
C     SAVE NVEC AS NOVEC
C
      NOVEC=NVEC
      ICK=0
C
C     SCHMIDT ORTHONORMALIZATION
C
      DO 100 K=1,NDSR*2
C
         DO ISTAT=1,NVEC
            BQA=DDOT(L7A,BVA(1,ISTAT),1,QA(1,K),1)
            BQB=DDOT(L7B,BVB(1,ISTAT),1,QB(1,K),1)
            BQ = BQA+BQB
            CALL DAXPY(L7A,-BQ,BVA(1,ISTAT),1,QA(1,K),1)
            CALL DAXPY(L7B,-BQ,BVB(1,ISTAT),1,QB(1,K),1)
         END DO
C
         FNORA =DDOT(L7A,QA(1,K),1,QA(1,K),1)
         FNORB =DDOT(L7B,QB(1,K),1,QB(1,K),1)
         FNOR  =FNORA+FNORB
         FNOR  =SQRT(FNOR)
         IF(FNOR.LT.THREB) GO TO 100
         IF(DBUG) THEN
            WRITE(6,'(A,I5,3F12.6)') 'FNORM', K, FNOR, FNORA, FNORB
         END IF
C
C     --- APPEND NEW B VECTOR
C
         NVEC=NVEC+1
         IF(K.LE.NDSR) THEN
            IINDX(NVEC)=IINDX(K)
         ELSE
            IINDX(NVEC)=IINDX(K-NDSR)
         END IF
         IF(NVEC.GT.MXVEC) THEN
            ICK=1
            RETURN
         END IF
         CALL VCLR(BVA(1,NVEC),1,L7)
         CALL VCLR(BVB(1,NVEC),1,L7)
         CALL DSCAL(L7A,1.0D+00/FNOR,QA(1,K),1)
         CALL DCOPY(L7A,QA(1,K),1,BVA(1,NVEC),1)
         CALL DSCAL(L7B,1.0D+00/FNOR,QB(1,K),1)
         CALL DCOPY(L7B,QB(1,K),1,BVB(1,NVEC),1)
C
 100  CONTINUE
C
C     ERROR TERMINATION
C
      MS=NVEC-NOVEC
      IF(MS.EQ.0) THEN
         IF(MASWRK) WRITE(IW,*)
     *        'SOMETHING WRONG...THERE ARE NO VECTORS TO BE ADDED'
         CALL ABRT
         STOP
      ENDIF
C
      IF(DBUG) THEN
      WRITE(IW,*)'NUMBER OF VECTORS ',NVEC,'=',NOVEC,'(OLD)+',MS,'(NEW)'
      ENDIF
      RETURN
      END
C
C*MODULE UTDDFT  *DECK URPAREDMS
      SUBROUTINE URPAREDMS(BVA,BVB,PMOA,MMOA,PMOB,MMOB,
     *                     XM_P,XM_M,NVEC,NOCA,NOCB,LX,L7)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DOUBLE PRECISION MMOA, MMOB
C
      DIMENSION BVA(L7,*)
      DIMENSION BVB(L7,*)
      DIMENSION PMOA(L7,*)
      DIMENSION MMOA(L7,*)
      DIMENSION PMOB(L7,*)
      DIMENSION MMOB(L7,*)
      DIMENSION XM_P(NVEC,*)
      DIMENSION XM_M(NVEC,*)
C
C     CONSTRUCT SMALL REDUCED (A+B) AND (A-B) MATRICES
C
      L7A=NOCA*(LX-NOCA)
      L7B=NOCB*(LX-NOCB)
      DO IP=1,NVEC
         DO IQ=1,NVEC
C
            DUMA=DDOT(L7A,BVA(1,IP),1,PMOA(1,IQ),1)
            DUMB=DDOT(L7B,BVB(1,IP),1,PMOB(1,IQ),1)
            XM_P(IP,IQ)=DUMA+DUMB
C
            DUMA=DDOT(L7A,BVA(1,IP),1,MMOA(1,IQ),1)
            DUMB=DDOT(L7B,BVB(1,IP),1,MMOB(1,IQ),1)
            XM_M(IP,IQ)=DUMA+DUMB
C
         ENDDO
      END DO
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK URPAEXPNDV
      SUBROUTINE URPAEXPNDV(VR,VL,VRO,VLO,BVA,BVB,
     *                      NDSR,NVEC,L7A,L7B,L7)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION VR(NVEC,*)
      DIMENSION VL(NVEC,*)
      DIMENSION VRO(L7,NDSR,2)
      DIMENSION VLO(L7,NDSR,2)
      DIMENSION BVA(L7,*),BVB(L7,*)
C
C     EXPAND REDUCED VECTORS TO REAL SIZE SPACE
C     VR AND VL : EIGENVECTOR
C     BVA, BVB  : TRIAL VECTOR
C
      CALL VCLR(VRO,1,L7*NDSR*2)
      CALL VCLR(VLO,1,L7*NDSR*2)
C
      DO N=1,NDSR
         DO JVEC=1,NVEC
            CALL DAXPY(L7A,VR(JVEC,N),BVA(1,JVEC),1,VRO(1,N,1),1)
            CALL DAXPY(L7A,VL(JVEC,N),BVA(1,JVEC),1,VLO(1,N,1),1)
            CALL DAXPY(L7B,VR(JVEC,N),BVB(1,JVEC),1,VRO(1,N,2),1)
            CALL DAXPY(L7B,VL(JVEC,N),BVB(1,JVEC),1,VLO(1,N,2),1)
         END DO
      END DO
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK URPAEXPNDV2
      SUBROUTINE URPAEXPNDV2(VR,VL,W_R,W_L,PMOA,PMOB,MMOA,MMOB,
     *                       NDSR,NVEC,L7A,L7B,L7)
C
C     EXPAND REDUCED VECTORS TO REAL SIZE SPACE
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DOUBLE PRECISION MMOA, MMOB
C
      DIMENSION VR(NVEC,*)
      DIMENSION VL(NVEC,*)
      DIMENSION W_R(L7,NDSR,2)
      DIMENSION W_L(L7,NDSR,2)
      DIMENSION PMOA(L7,*),PMOB(L7,*)
      DIMENSION MMOA(L7,*),MMOB(L7,*)
C
C     WA_R= PMOAA * VR
C     WB_R= PMOBB * VR
C     WA_L= MMOAA * VL
C     WB_L= MMOBB * VL
C
      CALL VCLR(W_R,1,L7*NDSR*2)
      CALL VCLR(W_L,1,L7*NDSR*2)
      DO N=1,NDSR
         DO JVEC=1,NVEC
            CALL DAXPY(L7A,VR(JVEC,N),PMOA(1,JVEC),1,W_R(1,N,1),1)
            CALL DAXPY(L7A,VL(JVEC,N),MMOA(1,JVEC),1,W_L(1,N,1),1)
            CALL DAXPY(L7B,VR(JVEC,N),PMOB(1,JVEC),1,W_R(1,N,2),1)
            CALL DAXPY(L7B,VL(JVEC,N),MMOB(1,JVEC),1,W_L(1,N,2),1)
         END DO
      END DO
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK URPAECHK
      SUBROUTINE URPAECHK(EE,INDX,IINDX,NVEC,NDSR,IMAX)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION EE(NVEC),INDX(NVEC),IINDX(NVEC)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      PARAMETER (TOEV=27.21138386D+00)
      PARAMETER (ZERO=0.0D+00)
C
C     REMOVE MINUS EIGENVALUES AND PRINT CURRENT EXCITATION ENERGIES
C
C     --- DEFINE INDX(NDIM)
C
      DO I=1,NVEC
         INDX(I)=I
      ENDDO
C
C     --- NUMBER OF MINUS EIGENVALUES : IMAX
C
      IMAX=0
      DO I=1,NVEC
         IF(EE(I).LT.ZERO) IMAX=I
      ENDDO
      IF(IMAX.NE.0) WRITE(IW,*) 'NUMBER OF MINUS EIGENVALUES=',IMAX
C
C     -- PRINT CURRENT EXCITATION ENERGIES
C
      DO ISTAT=1,NDSR
         EE(ISTAT)=SQRT(ABS(EE(ISTAT)))
         YDUM=EE(ISTAT)*TOEV !(627.5095D0/23.06035D0)
         IF(MASWRK) THEN
            IF(IINDX(ISTAT).EQ.1)
     *           WRITE(IW,8888) ISTAT, YDUM, '(A)'
            IF(IINDX(ISTAT).EQ.2)
     *           WRITE(IW,8888) ISTAT, YDUM, '(B)'
         END IF
      ENDDO
C
      RETURN
 8888 FORMAT(1X,'   STATE NUMBER',I4,2X,'ENERGY =',F12.6,1X,'EV',A5)
      END
C
C*MODULE UTDDFT  *DECK UTDOSCALC
      SUBROUTINE UTDOSCALC(OS,TXYZ,TDENA,TDENB,SCR,VRO,VLO,EE,VA,VB,
     *     AX,AY,AZ,L1,LX,L2,L7,L7A,L7C,NDSR,INXOV)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION VRO(L7,NDSR,2)
      DIMENSION VLO(L7,NDSR,2)
      DIMENSION EE(*)
      DIMENSION VA(L1,LX), VB(L1,LX)
      DIMENSION TDENA(L1,L1),TDENB(L1,L1),SCR(L1,L1)
      DIMENSION AX(L2),AY(L2),AZ(L2)
      DIMENSION TXYZ(3,NDSR)
      DIMENSION OS(NDSR)
      DIMENSION INXOV(2,*) ! (2,L7C)
C
      PARAMETER (MXAO=8192)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /XYZPRP/ XP,YP,ZP,
     *                DMX,DMY,DMZ,
     *                QXX,QYY,QZZ,QXY,QXZ,QYZ,
     *                QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ,
     *                OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                OXZZ,OYZZ,OZZZ,OXYZ,
     *                OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY,
     *                OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ
C
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
C
      DATA      HALF /0.5D+00/
C
C     --- CALCULATE OSCILLATOR STRENGTH
C
C        GET DIPOLE INTEGRALS AT CENTER OF MASS
C
      CALL CALCOM(XP,YP,ZP)
      CALL DIPINT(XP,YP,ZP,.FALSE.)
      CALL DAREAD(IDAF,IODA,AX,L2,95,0)
      CALL DAREAD(IDAF,IODA,AY,L2,96,0)
      CALL DAREAD(IDAF,IODA,AZ,L2,97,0)
C
      DO IST=1,NDSR
C
C     -- TRANSITION DENSITY
C
         CALL UTRAD(TDENA,TDENB,VRO,VLO,LX,L7,L7A,L7C,INXOV,NDSR,IST)
C
C     -- TRANSFORM TO AO-BASIS
C     --- ALPHA ---
C     TDENA = VA * TDENA * VA^T
         CALL DGEMM('N','N',L1,LX,LX,ONE,VA,L1,TDENA,LX,ZERO,SCR,L1)
         CALL DGEMM('N','T',L1,L1,LX,ONE,SCR,L1,VA,L1,ZERO,TDENA,L1)
C     --- BETA ---
C     TDENB = VB * TDENB * VB^T
         CALL DGEMM('N','N',L1,LX,LX,ONE,VB,L1,TDENB,LX,ZERO,SCR,L1)
         CALL DGEMM('N','T',L1,L1,LX,ONE,SCR,L1,VB,L1,ZERO,TDENB,L1)
C
C     -- TRANSITION DIPOLE
         DUMNX=ZERO
         DUMNY=ZERO
         DUMNZ=ZERO
         DO M=1,L1
            DO N=1,L1
               MN=IA(MAX(M,N))+MIN(M,N)
               TDEN= TDENA(M,N)+TDENB(M,N)
               DUMX=TDEN*AX(MN)
               DUMY=TDEN*AY(MN)
               DUMZ=TDEN*AZ(MN)
               DUMNX=DUMNX+DUMX
               DUMNY=DUMNY+DUMY
               DUMNZ=DUMNZ+DUMZ
            END DO
         END DO
         TXYZ(1,IST)=-DUMNX*HALF
         TXYZ(2,IST)=-DUMNY*HALF
         TXYZ(3,IST)=-DUMNZ*HALF
C
      ENDDO
C
C     -- OSCILLATOR STRENGTH, AND PRINTING
C
      CALL OSCALC(OS,TXYZ,EE,NDSR)
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK UTDPRNT
      SUBROUTINE UTDPRNT(OS,LABSYM,TXYZ,INXOV,NDSR,EE,VRO,IRMOA,IRMOB,
     *     L7,L7A,L7C,L0,IMAX)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,TRIPLET,CHECK,ABEL,ABELPT,
     *        SKIP,SG1T
C
      PARAMETER (MXATM=2000, MXRT=100)
C
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FUNCT / ENERGY,EG(3*MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INFOTD/ TDM(3,MXRT),CNVTOL,TRIPLET,MOCC(MXRT),MVIR(MXRT),
     *                NRADT,NTHET,NPHIT,NSTAT,NTRIAL,MAXVEC,NTHST,
     *                IRECTD,NLEBT,SG1T
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      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
      DIMENSION VRO(L7,NDSR,2),INXOV(2,*) ! (2,L7C)
      DIMENSION OS(NDSR),LABSYM(NDSR),TXYZ(3,NDSR),EE(NDSR)
      DIMENSION IRMOA(L0),IRMOB(L0)
C
      PARAMETER (TOEV=27.21138386D+00, THRESHOLD=3.0D-02)
      DATA IQMARK,IBLANK/4H????,4H    /
C
      IF(MASWRK) WRITE(IW,8885) MUL
C
      ABEL = ABELPT()
C
C              MIKE SAYS:
C         THE APPROACH USED FOR THE GROUND STATE, NAMELY ASSUME THE
C         STATE SYMMETRY OF THE GROUND STATE IS TOTALLY SYMMETRIC,
C         AND THUS GET EXCITED STATE SYMMETRY BY TAKING JUST THE
C         DIRECT PRODUCT OF THE LOWER AND UPPER ORBITAL OF EACH
C         EXCITATION IS INCORRECT.
C         THEREFORE -SKIP- IS INTRODUCED IN ORDER TO PRINT BLANK
C         SPACES FOR EACH STATE.
C         ASSUMING THAT WERE FIXED, PERHAPS BY ACCUMULATING THE
C         GROUND STATE SYMMETRY INTO -KRP-, AND EXTENDING THE
C         DIRECT PRODUCT TO IRP TIMES JRP TIMES KRP, WE MIGHT
C         WANT TO GO ON TO...
C         SEE ALSO COMMENTS IN THE CLOSED SHELL ROUTINE ABOUT THE
C         POSSIBILITY OF DOING SOMETHING IN NON-ABELIAN GROUPS TOO.
C
      SKIP = .TRUE.
C
      DO ISTAT=1,NDSR
         IF(ISTAT.LE.IMAX) GOTO 100
         OSI=OS(ISTAT)
         YDUM = TOEV*EE(ISTAT)
C
C      -- SYMMETRY OF EXCITED STATE
C
         IF(ABEL) THEN
C              REMOVE NEXT THREE LINES, AND THE MATCHING ENDIF,
C              AND ALL OTHER OCCURENCES OF -SKIP-,
C              IF YOU ARE THE PERSON TRYING TO FIX THIS ROUTINE.
            IF(SKIP) THEN
               LABSYM(ISTAT)=IBLANK
            ELSE
            CHECK = .FALSE.
            NAMSAV = -179
            DO IJ=1,L7C
               IF (IJ.LE.L7A) THEN
                  I=INXOV(1,IJ)
                  J=INXOV(2,IJ)
                  XDUM=VRO(IJ,ISTAT,1)
                  IRP = IRMOA(I)
                  JRP = IRMOA(J)
               ELSE
                  II=IJ-L7A
                  I=INXOV(1,IJ)
                  J=INXOV(2,IJ)
                  XDUM=VRO(II,ISTAT,2)
                  IRP = IRMOB(I)
                  JRP = IRMOB(J)
               END IF
               IF (ABS(XDUM).GT.THRESHOLD) THEN
                  IINAM = IJREP(2,1,IRP,JRP)
                  IF(CHECK) THEN
                     IF(IINAM.NE.NAMSAV) THEN
                        IF(MASWRK) WRITE(IW,8892)
                        LABSYM(ISTAT) = IQMARK
                        GO TO 50
                     END IF
                  ELSE
                     CHECK = .TRUE.
                     LABSYM(ISTAT) = IRPNAM(IINAM)
                     NAMSAV = IINAM
                  END IF
               END IF
            ENDDO
            END IF
         ELSE
            LABSYM(ISTAT)=IBLANK
         END IF
C
 50      CONTINUE
         IF(MASWRK) THEN
            WRITE(IW,8888) ISTAT, YDUM
            WRITE(IW,8889) OSI
            IF(.NOT.SKIP) WRITE(IW,8890) LABSYM(ISTAT)
            WRITE(IW,8887)
            WRITE(IW,8884)
         ENDIF
C
C            PRINT AMPLITUDE INFORMATION
C
         DO IJ=1,L7C
            IF(IJ.LE.L7A) THEN
C     ALPHA EXCITATIONS
               I=INXOV(1,IJ)
               J=INXOV(2,IJ)
               XDUM=VRO(IJ,ISTAT,1)
               IF (ABS(XDUM).GT.THRESHOLD) THEN
                  IF(MASWRK) WRITE(IW,8891) '(A)',IJ,XDUM,I,J
               ENDIF
            ELSE
C     BETA EXCITATIONS
               II=IJ-L7A
               I=INXOV(1,IJ)
               J=INXOV(2,IJ)
               XDUM=VRO(II,ISTAT,2)
               IF (ABS(XDUM).GT.THRESHOLD) THEN
                  IF(MASWRK) WRITE(IW,8891) '(B)',IJ,XDUM,I,J
               ENDIF
            END IF
         ENDDO
 100     CONTINUE
      END DO
C
C         PRINT SUMMARY TABLE
C
      IF(MASWRK) THEN
         WRITE(IW,9010) MUL
         IF(SKIP) THEN
            WRITE(IW,9020) 0,IBLANK,ESCF,0.0D+00
         ELSE
         WRITE(IW,9020) 0,IRPNAM(1),ESCF,0.0D+00
         END IF
         DO IST=1,NDSR
            WRITE(IW,9020) IST,LABSYM(IST),ESCF+EE(IST),TOEV*EE(IST),
     *                     (TXYZ(II,IST),II=1,3),OS(IST)
         ENDDO
         WRITE(IW,*) ' '
      END IF
C
C                ----- SAVE ENERGY RESULTS -----
C     CAN'T UPDATE K.E., ET CETERA WITHOUT A DENSITY MATRIX
C
      ENERGY    = ESCF + EE(NTHST)
      ETOT      = ENERGY
C     ESTATE(1) = ENERGY
      STATN     = MIN(1+NDSR,MXRT)
      DO I=1,MIN(NDSR,MXRT-1)
        ESTATE(I+1) = ESCF + EE(NTHST)
      ENDDO
      NSTMAX=NDSR
      DO I=1,NSTMAX
      ESTATE(I)       = EE(I)
      TDM(1,I)    = TXYZ(1,I)
      TDM(2,I)    = TXYZ(2,I)
      TDM(3,I)    = TXYZ(3,I)
      ENDDO
      CALL DAWRIT(IDAF,IODA,ENUCR,MXRT+15,2,0)
      IF(MASWRK) WRITE(IW,9030) NTHST,ENERGY
C
      RETURN
C
 8892 FORMAT(/1X,'WARNING: THIS STATE HAS BROKEN SYMMETRY, CHECK MOS')
 8890 FORMAT(1X,'SYMMETRY OF STATE   =',4X,A4)
 8889 FORMAT(1X,'OSCILLATOR STRENGTH =',F12.6)
 8888 FORMAT(/1X,'STATE #',I4,2X,'ENERGY =',F12.6,1X,'EV')
 8887 FORMAT(8X,'SPIN',4X,'DRF',4X,'COEF',7X,'OCC',6X,'VIR')
 8885 FORMAT(/10X,50(1H-)/
     *        10X,'RESULTS FOR EXCITED STATES WITH SPIN MULTIPLITY',I3/
     *        10X,50(1H-))
 8884 FORMAT(8X,4(1H-),4X,3(1H-),4X,4(1H-),7X,3(1H-),6X,3(1H-))
 8891 FORMAT(7X,A4,3X,I5,1X,F9.6,4X,I4,2X,'->',1X,I4,2X)
 9010 FORMAT(/18X,'SUMMARY OF TDDFT RESULTS, MULTIPLICITY=',I3//
     *         3X,'STATE',13X,'ENERGY',5X,'EXCITATION',
     *         6X,'TRANSITION DIPOLE',4X,'OSCILLATOR'/
     *        20X,'HARTREE',10X,'EV',9X,'X',7X,'Y',7X,'Z',5X,'STRENGTH')
 9020 FORMAT(1X,I3,2X,A4,1X,F20.10,1X,F8.3,3X,3F8.4,1X,F8.3)
 9030 FORMAT(1X,'SELECTING EXCITED STATE IROOT=',I4,' AT E=',F20.10/
     *       1X,'AS THE STATE OF INTEREST.'/)
      END
C
C*MODULE UTDDFT  *DECK UTRAD
      SUBROUTINE UTRAD(TDENA,TDENB,VRO,VLO,LX,L7,L7A,L7C,INXOV,NDSR,IST)
C
C     -- TRANSITION DENSITY
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION TDENA(LX,LX),TDENB(LX,LX)
      DIMENSION VRO(L7,NDSR,2)
      DIMENSION VLO(L7,NDSR,2)
      DIMENSION INXOV(2,*)
C     --CLEAR
      CALL VCLR(TDENA,1,2*LX*LX)
      CALL VCLR(TDENB,1,2*LX*LX)
C
      DO IPR=1,L7C
         I=INXOV(1,IPR)
         J=INXOV(2,IPR)
         IF(IPR.LE.L7A) THEN
C     --- ALPHA ORBITAL ---
C     -- IA COMPONENTS
            TDENA(I,J)=VRO(IPR,IST,1)+VLO(IPR,IST,1)
C     -- AI COMPONENTS
            TDENA(J,I)=VRO(IPR,IST,1)-VLO(IPR,IST,1)
         ELSE
            IJ = IPR-L7A
C     -- IA COMPONENTS
            TDENB(I,J)=VRO(IJ,IST,2)+VLO(IJ,IST,2)
C     -- AI COMPONENTS
            TDENB(J,I)=VRO(IJ,IST,2)-VLO(IJ,IST,2)
         END IF
      ENDDO
C
      RETURN
      END
C
C  HERE DOWN IS FUNCTIONAL SUPPORT
C
C*MODULE UTDDFT  *DECK UTDDFTSET
      SUBROUTINE UTDDFTSET(XYZGRD,XYZWGT,XYZDCH,VA,VB,
     *                     DRHOI_,AOMAX,GMO,ILENG,NPTGRD,L1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL DBUG
      DIMENSION VA(*),VB(*)
      DIMENSION DRHOI_(*)
      DIMENSION XYZGRD(3,NPTGRD),XYZWGT(NPTGRD)
      DIMENSION XYZDCH(*)
      DIMENSION AOMAX(L1)
      DIMENSION GMO(*)
      DBUG =.FALSE.
C
C     ---- READ GRID INFORMATION
C
      CALL GRDRD(22,XYZGRD,XYZWGT,IDUMMY,NPTGRD,0)
      IF (DBUG) THEN
      WRITE(6,*)'IN TDDFTSET'
      CALL TIMIT(1)
      ENDIF
C
C     ---- DENSITY
C
      CALL UDENMEM(RCUTOFF,RHO,DRHOI_,VA,VB,
     *     XYZWGT,XYZGRD,XYZDCH,AOMAX,GMO,NPTGRD,L1,ILENG)
      IF (DBUG) THEN
         WRITE(6,*)'DENMEM DONE'
         WRITE(6,*)'INTEGRATED RHO=',RHO
         CALL TIMIT(1)
      ENDIF
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK UDENMEM
      SUBROUTINE UDENMEM(RCUTOFF,RHO,DRHOI_,VA,VB,XYZWGT,XYZGRD,
     *                   XYZDCH,AOMAX,GMO,NPTGRD,L1,ILENG)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DRHOI_(*),VA(L1,*),VB(L1,*),
     *          XYZWGT(*),XYZGRD(3,*),XYZDCH(4,L1,*),AOMAX(*),GMO(*)
C
      LOGICAL SG1,GOPARR,DSKWRK,MASWRK,NXT
C
      PARAMETER (MXATM=2000, MXGRID=10, MXGRIDTYP=10)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /DFPRUN/ PRUNERADS(MXGRID,MXGRIDTYP),
     *                PRUNEATOMS(2,MXGRIDTYP),
     *                IPRUNECUTS(MXATM),NTOTGRIDPOINTS(MXATM),
     *                NGRIDS,MAXANG,NGRIDTYPS
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
      DATA ZERO      /0.0D+00/
C
C     SET CUT-OFFS
C
C      NPT=NRAD*NTHE*NPHI
      NPT=NRAD*MAXANG
      RCUTOFF=CONVHF/(NPT*NAT)
      CCUTOFF=1.0D-03/(NPT*NAT)
      IF (DFTTHR.EQ.ZERO) DFTTHR=1.0D-04/(NPT*NAT)
C
      IST=0
      IEND=0
      IAREA=0
C
      CALL VCLR(DRHOI_,1,2*4*NPTGRD)  ! FOR (A)&(B) SPIN
C
      RHO=ZERO
C
C     LOAD BALANCING TERMS
C
      LOOP = 0
      NGRAN = 20
      NLENG = (NPTGRD - 1)/ILENG + 1
      MCHUNK = (NLENG - 1)/(NGRAN*NPROC) + 1
      NEXT   = -1
      NXT    = IBTYP.EQ.1
C
C     CHANGED GOTO INTO A LOOP FOR PARALLELIZATION
C
      DO 10 ILG=1,NLENG
C
      IF(GOPARR) THEN
         IF(NXT) THEN
            LOOP = LOOP + 1
            ICHUNK=(LOOP-1)/MCHUNK
            IF(ICHUNK.GT.NEXT) THEN
               CALL DDI_DLBNEXT(NEXT)
            ENDIF
            IF(NEXT.NE.ICHUNK) GOTO 10
         ELSE
            IF(MOD(ILG,NPROC).NE.ME) GOTO 10
         ENDIF
      ENDIF
      IST =(ILG - 1)*ILENG + 1
      IEND=IST+ILENG-1
      IAREA=IAREA+1
      IF(IEND.GT.NPTGRD) IEND=NPTGRD
C************************************************************
C
C     -- MAKE AO-GRID & AOMAX
      CALL GRDGAO(IST,IEND,XYZGRD,XYZDCH,AOMAX,L1)
C
C     -- NUMERICAL INTEGRATION
C
      CALL UDENCNST(RCUTOFF,CCUTOFF,RHO,DRHOI_,VA,VB,XYZWGT,XYZDCH,
     *     GMO,IST,IEND,L1,NA,NB,NPTGRD)
C************************************************************
C
 10   CONTINUE
C
      IF (GOPARR) THEN
         IF (NXT) CALL DDI_DLBRESET
         CALL DDI_GSUMF(2350,DRHOI_,2*4*NPTGRD)
         CALL DDI_GSUMF(2351,RHO,1)
         IF(MASWRK) PRINT *, "RHO  ", RHO
      END IF
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK UDENCNST
      SUBROUTINE UDENCNST(RCUTOFF,CCUTOFF,RHO,DRHOI_,VA,VB,WGT,GAO,
     *                    GMO,IST,IEND,L1,NA,NB,NPTGRD)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION DRHOI_(NPTGRD,4,2)
      DIMENSION VA(L1,*), VB(L1,*)
      DIMENSION WGT(NPTGRD)
      DIMENSION GAO(4,L1,*)
      DIMENSION GMO(L1,4,2)
      DATA TWO       /2.0D+00/
C
C     -- GRID LOOP
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
C
C     -- TRANSFORM
C
         CALL VCLR(GMO,1,8*L1) ! FOR (A) AND (B) SPINS
         DO I=1,L1
            AO=GAO(1,I,IIPT)
            IF(ABS(AO).GT.CCUTOFF) THEN
               GX=GAO(2,I,IIPT)
               GY=GAO(3,I,IIPT)
               GZ=GAO(4,I,IIPT)
C     FOR (A) SPIN
               DO J=1,NA
                  COA=VA(I,J)
                  GMO(J,1,1)=GMO(J,1,1)+AO*COA
                  GMO(J,2,1)=GMO(J,2,1)+GX*COA
                  GMO(J,3,1)=GMO(J,3,1)+GY*COA
                  GMO(J,4,1)=GMO(J,4,1)+GZ*COA
               END DO
C     FOR (B) SPIN
               DO J=1,NB
                  COB=VB(I,J)
                  GMO(J,1,2)=GMO(J,1,2)+AO*COB
                  GMO(J,2,2)=GMO(J,2,2)+GX*COB
                  GMO(J,3,2)=GMO(J,3,2)+GY*COB
                  GMO(J,4,2)=GMO(J,4,2)+GZ*COB
               END DO
            ENDIF
         END DO
C
C     -- FORM DENSITY AT THIS POINT
C
         ROA=DDOT(NA,GMO(1,1,1),1,GMO(1,1,1),1)
         ROB=DDOT(NB,GMO(1,1,2),1,GMO(1,1,2),1)
C         ROB=ROA
         IF(ROA+ROB.LT.RCUTOFF) THEN
            GOTO 10
         ENDIF
C     FOR (A) SPIN
         GRADXA=TWO*DDOT(NA,GMO(1,1,1),1,GMO(1,2,1),1)
         GRADYA=TWO*DDOT(NA,GMO(1,1,1),1,GMO(1,3,1),1)
         GRADZA=TWO*DDOT(NA,GMO(1,1,1),1,GMO(1,4,1),1)
C
         DRHOI_(IPT,1,1)=ROA
         DRHOI_(IPT,2,1)=GRADXA
         DRHOI_(IPT,3,1)=GRADYA
         DRHOI_(IPT,4,1)=GRADZA
C     FOR (B) SPIN
         GRADXB=TWO*DDOT(NB,GMO(1,1,2),1,GMO(1,2,2),1)
         GRADYB=TWO*DDOT(NB,GMO(1,1,2),1,GMO(1,3,2),1)
         GRADZB=TWO*DDOT(NB,GMO(1,1,2),1,GMO(1,4,2),1)
C
         DRHOI_(IPT,1,2)=ROB
         DRHOI_(IPT,2,2)=GRADXB
         DRHOI_(IPT,3,2)=GRADYB
         DRHOI_(IPT,4,2)=GRADZB
C
C         RHO=RHO+TWO*ROA*WGT(IPT)
         RHO = RHO + (ROA+ROB)*WGT(IPT)
 10      CONTINUE
      END DO
C     -- GRID LOOP END
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK UTDFXCP
      SUBROUTINE UTDFXCP(FXC,RHO,XYZGRD,XYZWGT,XYZDCH,
     *     VHA,VHB,VA,VB,DRHOI_,DTRAI_,
     *     COEF,EX,EC,EX0,EC0,AOMAX,GMO,VALGA,ILENG,NPTGRD,
     *     L1,L2,LX,NDER)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION FXC(L2,2)
      DIMENSION VHA(LX,L1),VHB(LX,L1)
      DIMENSION VA(L1,*),VB(L1,*)
      DIMENSION XYZGRD(3,*)
      DIMENSION XYZWGT(*)
      DIMENSION XYZDCH(*)
      DIMENSION DRHOI_(NPTGRD,4,2)
      DIMENSION DTRAI_(8,*)
      DIMENSION COEF  (4,4,*) ! 1=>AA, 2=>AB, 3=>BA, 4=>BB
      DIMENSION EX(*), EC(*), EX0(*), EC0(*)
      DIMENSION AOMAX(L1)
      DIMENSION GMO(*)          ! GMO(L1,4,2)
      DIMENSION VALGA(*)        ! VALGA(L1,2)
C
      PARAMETER (MXATM=2000, MXGRID=10, MXGRIDTYP=10)
C
      LOGICAL DBUG,SG1,NXT,GOPARR,DSKWRK,MASWRK
C
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /DFPRUN/ PRUNERADS(MXGRID,MXGRIDTYP),
     *                PRUNEATOMS(2,MXGRIDTYP),
     *                IPRUNECUTS(MXATM),NTOTGRIDPOINTS(MXATM),
     *                NGRIDS,MAXANG,NGRIDTYPS
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA ZERO      /0.0D+00/
C
      DBUG =.FALSE.
C     SET CUT-OFFS
C
C      NPT=NRAD*NTHE*NPHI
      NPT=NRAD*MAXANG
      RCUTOFF=CONVHF/(NPT*NAT)
      IF(NDER.EQ.3) RCUTOFF=1.0D-12
      CCUTOFF=1.0D-03/(NPT*NAT)
      IF (DFTTHR.EQ.ZERO) DFTTHR=1.0D-04/(NPT*NAT)
      IF(DBUG) THEN
         WRITE(6,*)'RCUTOFF/CCUTOFF=',RCUTOFF,CCUTOFF
      ENDIF
C
      IST=0
      IEND=0
      IAREA=0
C     
      CALL VCLR(FXC,1,L2*2)
      RHO=ZERO
C     FOR PARALLEL
      NEXT = -1
      MINE = -1
      IPCOUNT = ME-1
      NXT = IBTYP.EQ.1
C
   10 CONTINUE
      IST=IEND+1
      IEND=IST+ILENG-1
      IAREA=IAREA+1
      IF(IEND.GT.NPTGRD) IEND=NPTGRD
C     PARALL
      IF(GOPARR) THEN
         IF(NXT) THEN
            MINE=MINE+1
            IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
            IF(NEXT.NE.MINE) GOTO 110
         ELSE
            IPCOUNT=IPCOUNT+1
            IF(MOD(IPCOUNT,NPROC).NE.0) GOTO 110
         END IF
      END IF
C
C     ---- FUNCTIONAL ENERGY AND DERIVATIVES
C          EX(ILENG,10) & EC(ILENG,20)
      IRU=2
      CALL DFXCSET(EX0,EX,EC0,EC,DRHOI_,XYZWGT,RCUTOFF,EXC,ECOR,
     *     ILENG,NPTGRD,IRU,NDER,IST,IEND)
C
C     -- MAKE AO-GRID & AOMAX
C
      CALL GRDGAO(IST,IEND,XYZGRD,XYZDCH,AOMAX,L1)
C
C     VH = MATRIX OF PSI_I PSI_A : FOR ALPHA ORBITAL
C     GET DTRAI FOR ALPHA ORBITAL
      CALL UTRADEN(RCUTOFF,CCUTOFF,DTRAI_,DRHOI_,VA,VB,VHA,VHB,
     *     XYZDCH,GMO,IST,IEND,L1,LX,NA,NB,NPTGRD,ILENG)
C
C     ---- COEF(NPTGRD,4,4)
C
      CALL UTDCOEF(COEF,EX,EC,DRHOI_,DTRAI_,RCUTOFF,NPTGRD,
     *     IST,IEND,ILENG)
C
C     -- FXC(L2,2)
C
      CALL UNUMFK(RHO,RCUTOFF,DFTTHR,FXC,DRHOI_,DTRAI_,COEF,XYZWGT,
     *     XYZDCH,VALGA,IST,IEND,L1,L2,NPTGRD)
C      CALL UNUMINTF(RHO,FXC,XYZWGT,XYZGRD,XYZDCH,DRHOI_,
C     *     DTRAI_,COEF,AOMAX,VALGA,NPTGRD,L1,L2,ILENG)
C
 110  CONTINUE
      IF(IEND.EQ.NPTGRD) THEN
         GOTO 20
      ELSE
         GOTO 10
      ENDIF
   20 CONTINUE
C
      IF (GOPARR) THEN
         IF(NXT) CALL DDI_DLBRESET
         CALL DDI_GSUMF(2330,FXC,L2*2)
         CALL DDI_GSUMF(2331,RHO,1)
      END IF
C
      IF (DBUG.AND.MASWRK) THEN
      WRITE(6,*)'TDFXCP DONE'
      WRITE(6,*)'INTEGRATED RHO=',RHO
      CALL TIMIT(1)
      ENDIF
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK UTDCOEF
      SUBROUTINE UTDCOEF(COEF,EX,EC,DRHOI_,DTRAI_,RCUTOFF,NPTGRD,
     *     IST,IEND,ILENG)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /XDERIX/ NXDIM(3),KRA,KRB,KGA,KGB,KRARA,KRAGA,KGAGA,
     *                KRBRB,KRBGB,KGBGB,KRARARA,KRARAGA,KRAGAGA,KGAGAGA,
     *                KRBRBRB,KRBRBGB,KRBGBGB,KGBGBGB
      COMMON/ CDERIX/ NCDIM(3),IRA,IRB,IGA,IGB,IGC,IRARA,IRARB,
     *                IRAGA,IRAGB,IRAGC,IRBRB,IRBGA,IRBGB,IRBGC,IGAGA,
     *                IGAGB,IGAGC,IGBGB,IGBGC,IGCGC,
     *                IRARARA,IRARARB,IRARAGA,IRARAGB,IRARAGC,IRARBRB,
     *                IRARBGA,IRARBGB,IRARBGC,IRAGAGA,IRAGAGB,IRAGAGC,
     *                IRAGBGB,IRAGBGC,IRAGCGC,IRBRBRB,IRBRBGA,IRBRBGB,
     *                IRBRBGC,IRBGAGA,IRBGAGB,IRBGAGC,IRBGBGB,IRBGBGC,
     *                IRBGCGC,IGAGAGA,IGAGAGB,IGAGAGC,IGAGBGB,IGAGBGC,
     *                IGAGCGC,IGBGBGB,IGBGBGC,IGBGCGC,IGCGCGC
      DIMENSION COEF(4,4,*)
      DIMENSION DRHOI_(NPTGRD,4,2)
      DIMENSION DTRAI_(4,2,*)
      DIMENSION EX(ILENG,*) ! (ILENG,10)
      DIMENSION EC(ILENG,*) ! (ILENG,20)
      DATA TWO       /2.0D+00/
      DATA FOUR      /4.0D+00/
C
      CALL VCLR(COEF,1,ILENG*16)
C
      DO IPT=IST,IEND
         IIPT = IPT-IST+1
C     --- FOR ALPHA ORBITAL
         RHOI=DRHOI_(IPT,1,1)
         IF(RHOI.GT.RCUTOFF) THEN
            DRHOIX = DRHOI_(IPT,2,1) ! D RHO_ALPHA
            DRHOIY = DRHOI_(IPT,3,1)
            DRHOIZ = DRHOI_(IPT,4,1)
            DRHOJX = DRHOI_(IPT,2,2) ! D RHO_BETA
            DRHOJY = DRHOI_(IPT,3,2)
            DRHOJZ = DRHOI_(IPT,4,2)
            DUMNB  = DTRAI_(1,1,IIPT) ! \PSI_IA
            DXB    = DTRAI_(2,1,IIPT)
            DYB    = DTRAI_(3,1,IIPT)
            DZB    = DTRAI_(4,1,IIPT)
C
C     -- DUMGRB: \DEL \PSI_ALPHA * \DEL RHO_ALPHA
            DUMGRB=DXB*DRHOIX+DYB*DRHOIY+DZB*DRHOIZ
            DUMGRC=DXB*DRHOJX+DYB*DRHOJY+DZB*DRHOJZ
C
C     ---- CONTRACTION
C     -- EXCHANGE
C     -- FOR DUMN (PM * PN):OK
            XDUM1B=EX(IIPT,KRARA)*DUMNB+TWO*EX(IIPT,KRAGA)*DUMGRB
C     -- FOR DUMGR(DEL R * DEL PM*PN):OK
            XDUM2B=TWO*EX(IIPT,KRAGA)*DUMNB+FOUR*EX(IIPT,KGAGA)*DUMGRB
C     -- FOR DUMXP:OK
            XDUM3B=TWO*EX(IIPT,KGA)
C     -- CORRELATION
C     -- FOR DUMN (PM*PN):OK
            CDUM1BAA=EC(IIPT,IRARA)*DUMNB+TWO*EC(IIPT,IRAGA)*DUMGRB
     *           +EC(IIPT,IRAGC)*DUMGRC
C     -- FOR DUMGR(DEL RHO_A * DEL PM*PN):OKX
            CDUM2BAA=TWO*EC(IIPT,IRAGA)*DUMNB
     *           +FOUR*EC(IIPT,IGAGA)*DUMGRB
     *           +TWO*EC(IIPT,IGAGC)*DUMGRC
            CDUM2BAC=EC(IIPT,IRAGC)*DUMNB+TWO*EC(IIPT,IGAGC)*DUMGRB
     *           +EC(IIPT,IGCGC)*DUMGRC
C     -- FOR DUMXB:OK
            CDUM3BAA=TWO*EC(IIPT,IGA)
C     -- TOTAL
            COEF(1,1,IIPT)=XDUM1B+CDUM1BAA
            COEF(2,1,IIPT)=XDUM2B+CDUM2BAA ! DEL RHO_A * DEL (PM * PN)
            COEF(3,1,IIPT)=       CDUM2BAC ! DEL RHO_B * DEL (PM * PN)
            COEF(4,1,IIPT)=XDUM3B+CDUM3BAA
C     ============================
C     ==== FOR BETA(IA) ALPHA (JB)
C     ============================
C
C     FOR DUMN(PM*PN)
            CDUM1BAB=EC(IIPT,IRARB)*DUMNB+TWO*EC(IIPT,IRBGA)*DUMGRB
     *           +EC(IIPT,IRBGC)*DUMGRC
C     FOR DEL RHO_B * DEL (PM*PN)
            CDUM2BAB=TWO*EC(IIPT,IRAGB)*DUMNB
     *           +FOUR*EC(IIPT,IGAGB)*DUMGRB
     *           +TWO*EC(IIPT,IGBGC)*DUMGRC
C     FOR DEL RHO_A * DEL (PM*PN)
            CDUM2BAC=EC(IIPT,IRAGC)*DUMNB+TWO*EC(IIPT,IGAGC)*DUMGRB
     *           +EC(IIPT,IGCGC)*DUMGRC
C     --
            CDUM3BAB = EC(IIPT,IGC)
C     -- TOTAL
            COEF(1,3,IIPT) = CDUM1BAB
            COEF(2,3,IIPT) = CDUM2BAB
            COEF(3,3,IIPT) = CDUM2BAC  !<=== FOR BETA
            COEF(4,3,IIPT) = CDUM3BAB
         ENDIF
C
C
C     --- FOR BETA ORBITAL
         RHOI=DRHOI_(IIPT,1,2)
         IF(RHOI.GT.RCUTOFF) THEN
            DRHOIX=DRHOI_(IPT,2,2) ! D RHO_BETA
            DRHOIY=DRHOI_(IPT,3,2)
            DRHOIZ=DRHOI_(IPT,4,2)
            DRHOJX=DRHOI_(IPT,2,1) ! D RHO_ALPHA
            DRHOJY=DRHOI_(IPT,3,1)
            DRHOJZ=DRHOI_(IPT,4,1)
            DUMNB =DTRAI_(1,2,IIPT) ! \PSI_IA
            DXB   =DTRAI_(2,2,IIPT)
            DYB   =DTRAI_(3,2,IIPT)
            DZB   =DTRAI_(4,2,IIPT)
C
C     -- DUMGRB: \DEL \PSI_BETA * \DEL RHO_BETA
            DUMGRB=DXB*DRHOIX+DYB*DRHOIY+DZB*DRHOIZ
C     -- DUMGRC: \DEL \PSI_BETA * \DEL RHO_ALPHA
            DUMGRC=DXB*DRHOJX+DYB*DRHOJY+DZB*DRHOJZ
C
C     ---- CONTRACTION
C     ============================
C     ==== FOR BETA BETA
C     ============================
C     -- EXCHANGE
C     -- FOR DUMN (PM * PN):OK
            XDUM1B=EX(IIPT,KRBRB)*DUMNB+TWO*EX(IIPT,KRBGB)*DUMGRB
C     -- FOR DUMGR(DEL R * DEL PM*PN):OK
            XDUM2B=TWO*EX(IIPT,KRBGB)*DUMNB+FOUR*EX(IIPT,KGBGB)*DUMGRB
C     -- FOR DUMXP:OK
            XDUM3B=TWO*EX(IIPT,KGB)
C     -- CORRELATION
C     -- FOR DUMN (PM*PN):OK
            CDUM1BAA=EC(IIPT,IRBRB)*DUMNB+TWO*EC(IIPT,IRBGB)*DUMGRB
     *           +EC(IIPT,IRBGC)*DUMGRC
C     -- FOR DUMGR(DEL RHO_B * DEL PM*PN):OKX
            CDUM2BAA=TWO*EC(IIPT,IRBGB)*DUMNB
     *           +FOUR*EC(IIPT,IGBGB)*DUMGRB
     *           +TWO*EC(IIPT,IGBGC)*DUMGRC
C     -- FOR DUMGC(DEL RHO_A * DEL PM*PN):OKX
            CDUM2BAC=EC(IIPT,IRBGC)*DUMNB+TWO*EC(IIPT,IGBGC)*DUMGRB
     *           +EC(IIPT,IGCGC)*DUMGRC
C     -- FOR DUMXB:OK
            CDUM3BAA=TWO*EC(IIPT,IGB)
C     -- TOTAL:BB
            COEF(1,4,IIPT)=XDUM1B+CDUM1BAA
            COEF(2,4,IIPT)=XDUM2B+CDUM2BAA ! DEL RHO_B * DEL (PM * PN)
            COEF(3,4,IIPT)=       CDUM2BAC ! DEL RHO_A * DEL (PM * PN)
            COEF(4,4,IIPT)=XDUM3B+CDUM3BAA
C     ============================
C     ==== FOR ALPHA BETA(JB)
C     ============================
C
C
C     FOR DUMN(PM*PN)
            CDUM1BAB=EC(IIPT,IRARB)*DUMNB+TWO*EC(IIPT,IRAGB)*DUMGRB
     *           +EC(IIPT,IRAGC)*DUMGRC
C     FOR DEL RHO_A * DEL (PM*PN)
            CDUM2BAB=TWO*EC(IIPT,IRBGA)*DUMNB
     *           +FOUR*EC(IIPT,IGAGB)*DUMGRB
     *           +TWO*EC(IIPT,IGAGC)*DUMGRC
C     FOR DEL RHO_B * DEL (PM*PN)
            CDUM2BAC=EC(IIPT,IRBGC)*DUMNB+TWO*EC(IIPT,IGBGC)*DUMGRB
     *           +EC(IIPT,IGCGC)*DUMGRC
C     --
            CDUM3BAB = EC(IIPT,IGC)
C     -- TOTAL
            COEF(1,2,IIPT) = CDUM1BAB
            COEF(2,2,IIPT) = CDUM2BAB
            COEF(3,2,IIPT) = CDUM2BAC
            COEF(4,2,IIPT) = CDUM3BAB
         ENDIF
      END DO
      RETURN
      END
C
C*MODULE UTDDFT  *DECK UNUMFK
      SUBROUTINE UNUMFK(RHO,RCUTOFF,CUTOFF,FXC,DRHOI_,DTRAI_,COEF,WGT,
     *     GAO,VALGA,IST,IEND,L1,L2,NPTGRD)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION FXC   (L2,2) !
      DIMENSION DRHOI_(NPTGRD,4,2)
      DIMENSION DTRAI_(4,2,*)
      DIMENSION COEF  (4,4,*)
      DIMENSION WGT   (NPTGRD)
      DIMENSION GAO   (4,L1,*)
      DIMENSION VALGA (L1,2)
      DATA ZERO      /0.0D+00/
      DATA TWO       /2.0D+00/
      LOGICAL  LPRT
C
C     -- GRID LOOP
C
      LPRT=.FALSE.
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
         RHOI = DRHOI_(IPT,1,1) ! ALPHA
         RHOJ = DRHOI_(IPT,1,2) ! BETA
         IF((RHOI+RHOJ).LT.RCUTOFF) GOTO 200
         FTOTWT=WGT(IPT)
C     FOR AA
         C1AA=COEF(1,1,IIPT)
         C2AA=COEF(2,1,IIPT)
         C3AA=COEF(3,1,IIPT)
         C4AA=COEF(4,1,IIPT)
C     FOR AB
         C1AB=COEF(1,2,IIPT)
         C2AB=COEF(2,2,IIPT)
         C3AB=COEF(3,2,IIPT)
         C4AB=COEF(4,2,IIPT)
C
         DRHOIX=DRHOI_(IPT,2,1) ! \DEL \RHO FOR ALPHA ORBITAL
         DRHOIY=DRHOI_(IPT,3,1)
         DRHOIZ=DRHOI_(IPT,4,1)
         DRHOJX=DRHOI_(IPT,2,2) ! \DEL \RHO FOR ALPHA ORBITAL
         DRHOJY=DRHOI_(IPT,3,2)
         DRHOJZ=DRHOI_(IPT,4,2)
         DXB   =DTRAI_(2,1,IIPT) ! \DEL \PSI_I \PSI_A FOR ALPHA ORBITAL
         DYB   =DTRAI_(3,1,IIPT)
         DZB   =DTRAI_(4,1,IIPT)
         DXB2  =DTRAI_(2,2,IIPT) ! \DEL \PSI_I \PSI_A FOR ALPHA ORBITAL
         DYB2  =DTRAI_(3,2,IIPT)
         DZB2  =DTRAI_(4,2,IIPT)
C
         IND=1
         AOMAX=ZERO
         AMAXVAL=ZERO
         DO M=1,L1
            PM =GAO(1,M,IIPT)
            PMX=GAO(2,M,IIPT)
            PMY=GAO(3,M,IIPT)
            PMZ=GAO(4,M,IIPT)
C     FOR ALPHA ALPHA
            DUM1 = C1AA/TWO*PM
            DUM2 = C2AA*(DRHOIX*PMX+DRHOIY*PMY+DRHOIZ*PMZ)
            DUM3 = C3AA*(DRHOJX*PMX+DRHOJY*PMY+DRHOJZ*PMZ)
            DUM4 = C4AA*(DXB*PMX+DYB*PMY+DZB*PMZ)
            VALGA(M,1)=(DUM1+DUM2+DUM3+DUM4)*FTOTWT
C     FOR ALPHA BETA
            DUM1 = C1AB/TWO*PM
            DUM2 = C2AB*(DRHOIX*PMX+DRHOIY*PMY+DRHOIZ*PMZ)
            DUM3 = C3AB*(DRHOJX*PMX+DRHOJY*PMY+DRHOJZ*PMZ)
            DUM4 = C4AB*(DXB2*PMX+DYB2*PMY+DZB2*PMZ)
            !IF (DUM4 .NE. 0.0D0) PRINT *, DUM4
            VALGA(M,2)=(DUM1+DUM2+DUM3+DUM4)*FTOTWT
C
            VMA=VALGA(M,1)
            VMB=VALGA(M,2)
            VM =VMA+VMB
            IF(LPRT) PRINT *, 'FK ', VM, PM
            AOMAX=MAX(AOMAX,ABS(PM))
            AMAXVAL=MAX(AMAXVAL,ABS(VM))
            CVAL=ABS(PM*AMAXVAL)+ABS(VM*AOMAX)
            IF(CVAL.LT.CUTOFF) THEN
               IND=IND+M
               GOTO 100
            ENDIF
            DO N=1,M
               PN=GAO(1,N,IIPT)
               FXC(IND,1)=FXC(IND,1)+PM*VALGA(N,1)+VMA*PN  ! FOR FXC_AA
               FXC(IND,1)=FXC(IND,1)+PM*VALGA(N,2)+VMB*PN  ! FOR FXC_AB
               IND=IND+1
            END DO
 100        CONTINUE
         END DO
C
C     FOR BETA
C     BB
         C1AA=COEF(1,4,IIPT)
         C2AA=COEF(2,4,IIPT)
         C3AA=COEF(3,4,IIPT)
         C4AA=COEF(4,4,IIPT)
C     BA
         C1AB=COEF(1,3,IIPT)
         C2AB=COEF(2,3,IIPT)
         C3AB=COEF(3,3,IIPT)
         C4AB=COEF(4,3,IIPT)
C
         DRHOIX=DRHOI_(IPT,2,2) ! \DEL \RHO (BETA)
         DRHOIY=DRHOI_(IPT,3,2)
         DRHOIZ=DRHOI_(IPT,4,2)
         DRHOJX=DRHOI_(IPT,2,1) ! \DEL \RHO FOR ALPHA ORBITAL
         DRHOJY=DRHOI_(IPT,3,1)
         DRHOJZ=DRHOI_(IPT,4,1)
         DXB  =DTRAI_(2,2,IIPT)  ! \DEL \PSI_I \PSI_A FOR ALPHA ORBITAL
         DYB  =DTRAI_(3,2,IIPT)
         DZB  =DTRAI_(4,2,IIPT)
         DXB2 =DTRAI_(2,1,IIPT)  ! \DEL \PSI_I \PSI_A FOR ALPHA ORBITAL
         DYB2 =DTRAI_(3,1,IIPT)
         DZB2 =DTRAI_(4,1,IIPT)
C
         IND=1
         AOMAX=ZERO
         AMAXVAL=ZERO
         DO M=1,L1
            PM =GAO(1,M,IIPT)
            PMX=GAO(2,M,IIPT)
            PMY=GAO(3,M,IIPT)
            PMZ=GAO(4,M,IIPT)
C     FOR BETA BETA
            DUM1 = C1AA/TWO*PM
C              DEL RHO_B * DEL PMN
            DUM2 = C2AA*(DRHOIX*PMX+DRHOIY*PMY+DRHOIZ*PMZ)
C              DEL RHO_A * DEL PMN
            DUM3 = C3AA*(DRHOJX*PMX+DRHOJY*PMY+DRHOJZ*PMZ)
            DUM4 = C4AA*(DXB*PMX+DYB*PMY+DZB*PMZ)
            !IF (DUM4 .NE. 0.0D0) PRINT *, DUM4
            VALGA(M,1)=(DUM1+DUM2+DUM3+DUM4)*FTOTWT
C     FOR BETA ALPHA
            DUM1 = C1AB/TWO*PM
            DUM2 = C2AB*(DRHOIX*PMX+DRHOIY*PMY+DRHOIZ*PMZ)
            DUM3 = C3AB*(DRHOJX*PMX+DRHOJY*PMY+DRHOJZ*PMZ)
            DUM4 = C4AB*(DXB2*PMX+DYB2*PMY+DZB2*PMZ)
           ! IF (DUM4 .NE. 0.0D0) PRINT *, DUM4
            VALGA(M,2)=(DUM1+DUM2+DUM3+DUM4)*FTOTWT
C
            VMA=VALGA(M,1)
            VMB=VALGA(M,2)
            VM =VMA+VMB
            AOMAX=MAX(AOMAX,ABS(PM))
            AMAXVAL=MAX(AMAXVAL,ABS(VM))
            CVAL=ABS(PM*AMAXVAL)+ABS(VM*AOMAX)
            IF(CVAL.LT.CUTOFF) THEN
               IND=IND+M
               GOTO 120
            ENDIF
            DO N=1,M
               PN=GAO(1,N,IIPT)
               FXC(IND,2)=FXC(IND,2)+PM*VALGA(N,1)+VMA*PN ! FOR BB
               FXC(IND,2)=FXC(IND,2)+PM*VALGA(N,2)+VMB*PN ! FOR BA
               IND=IND+1
            END DO
 120        CONTINUE
         END DO
C
C     -- DEBUG DENSITY
         RHO   =RHO   + (RHOI+RHOJ) *WGT(IPT)
  200 CONTINUE
      END DO
C     --END GRID LOOP
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK UTRADEN
      SUBROUTINE UTRADEN(RCUTOFF,CCUTOFF,DTRAI_,DRHOI_,VA,VB,VHA,VHB,
     *                   GAO,GMO,IST,IEND,L1,LX,NA,NB,NPTGRD,ILENG)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION DRHOI_(NPTGRD,4,2)
C         DTRAI_(1:4,1,*) ARE ALPHA,   DTRAI_(1:4,2,*) ARE BETA
      DIMENSION DTRAI_(4,2,*)
      DIMENSION VA(L1,*),VB(L1,*)
      DIMENSION VHA(LX,*),VHB(LX,*)
      DIMENSION GAO(4,L1,*)
      DIMENSION GMO(L1,4,2)
      DATA ZERO      /0.0D+00/
C
      CALL VCLR(DTRAI_,1,8*ILENG)
C
C     -- GRID LOOP
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
         ROA=DRHOI_(IPT,1,1)
         ROB=DRHOI_(IPT,1,2)
         IF(ROA+ROB.LT.RCUTOFF) GOTO 10
C
C     -- TRANSFORM
C
         CALL VCLR(GMO,1,8*L1)
         DO I=1,L1
            AO=GAO(1,I,IIPT)
            IF(ABS(AO).GT.CCUTOFF) THEN
               GX=GAO(2,I,IIPT)
               GY=GAO(3,I,IIPT)
               GZ=GAO(4,I,IIPT)
C     FOR ALPHA ORBITAL
               DO J=1,NA
                  COA=VA(I,J)
                  GMO(J,1,1)=GMO(J,1,1)+AO*COA
                  GMO(J,2,1)=GMO(J,2,1)+GX*COA
                  GMO(J,3,1)=GMO(J,3,1)+GY*COA
                  GMO(J,4,1)=GMO(J,4,1)+GZ*COA
               END DO
C     FOR BETA ORBITAL
               DO J=1,NB
                  COB=VB(I,J)
                  GMO(J,1,2)=GMO(J,1,2)+AO*COB
                  GMO(J,2,2)=GMO(J,2,2)+GX*COB
                  GMO(J,3,2)=GMO(J,3,2)+GY*COB
                  GMO(J,4,2)=GMO(J,4,2)+GZ*COB
               END DO
            ENDIF
         END DO
C
C     -- TRANSITION DENSITY
C
         T1A=ZERO
         TXA=ZERO
         TYA=ZERO
         TZA=ZERO
         T1B=ZERO
         TXB=ZERO
         TYB=ZERO
         TZB=ZERO
         DO I=1,L1
            AOI=GAO(1,I,IIPT)
            IF(ABS(AOI).GT.CCUTOFF) THEN
               GXI=GAO(2,I,IIPT)
               GYI=GAO(3,I,IIPT)
               GZI=GAO(4,I,IIPT)
C
               D0A=ZERO
               DXA=ZERO
               DYA=ZERO
               DZA=ZERO
C     FOR ALPHA ORBITAL
               DO J=1,NA
                  VHAJI = VHA(J,I)
                  D0A=D0A + VHAJI*GMO(J,1,1)
                  DXA=DXA + VHAJI*GMO(J,2,1)
                  DYA=DYA + VHAJI*GMO(J,3,1)
                  DZA=DZA + VHAJI*GMO(J,4,1)
               END DO
C     FOR BETA ORBITAL
               D0B=ZERO
               DXB=ZERO
               DYB=ZERO
               DZB=ZERO
               DO J=1,NB
                  VHBJI = VHB(J,I)
                  D0B=D0B + VHBJI*GMO(J,1,2)
                  DXB=DXB + VHBJI*GMO(J,2,2)
                  DYB=DYB + VHBJI*GMO(J,3,2)
                  DZB=DZB + VHBJI*GMO(J,4,2)
               END DO
C
               T1A=T1A+D0A*AOI
               TXA=TXA+D0A*GXI+DXA*AOI
               TYA=TYA+D0A*GYI+DYA*AOI
               TZA=TZA+D0A*GZI+DZA*AOI
C
               T1B=T1B+D0B*AOI
               TXB=TXB+D0B*GXI+DXB*AOI
               TYB=TYB+D0B*GYI+DYB*AOI
               TZB=TZB+D0B*GZI+DZB*AOI
            ENDIF
         END DO
C
         DTRAI_(1,1,IIPT)=T1A
         DTRAI_(2,1,IIPT)=TXA
         DTRAI_(3,1,IIPT)=TYA
         DTRAI_(4,1,IIPT)=TZA
C
         DTRAI_(1,2,IIPT)=T1B
         DTRAI_(2,2,IIPT)=TXB
         DTRAI_(3,2,IIPT)=TYB
         DTRAI_(4,2,IIPT)=TZB
 10      CONTINUE
      END DO
C     -- GRID LOOP END
C
      RETURN
      END
C
C  HERE DOWN IS INTEGRAL SUPPORT
C
C*MODULE UTDDFT  *DECK UTD2E
      SUBROUTINE UTD2E(PA,PB,FA,FB,FM2A,FM2B,
     *                 BUF,IBUF,NBF,SCAL,IMF,IPF,NV)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,TDSKWRK
      LOGICAL PACK2E,SCHWRZ,DIRSCF,FDIFF,LCFLAG,LRINT
C
      DIMENSION PA(*),PB(*)
      DIMENSION FA(*),FB(*),FM2A(*),FM2B(*)
      DIMENSION BUF(*),IBUF(*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXAO=8192)
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      COMMON /TDSCL/ SCLINT,IMFINT,IPFINT
      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,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     &                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     &                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     &                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
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
      IMFINT = IMF
      IPFINT = IPF
      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
            CALL SHLTD(SCFTYP,PA,PB,XX(LDSH),NBF,NBF2,NV)
C           READ EXCHANGE-INTEGRAL SCREENING INFORMATION
            CALL DAREAD(IDAF,IODA,XX(LXINTS),NSH2,54,0)
         END IF
C
         IF(LCFLAG) THEN
         LRINT=.TRUE.
         CALL UTDTWOEI(SCFTYP,SCHWRZ,NINT,NSCHWZ,NBF,NBF2,
     *              XX(LXINTS),NSH2,XX(LGHOND),MAXG,
     *              IA,PA,PB,FA,FB,FM2A,FM2B,XX(LDSH),NV)
         LRINT=.FALSE.
         ELSE
         LRINT=.FALSE.
         END IF
         CALL UTDTWOEI(SCFTYP,SCHWRZ,NINT,NSCHWZ,NBF,NBF2,
     *              XX(LXINTS),NSH2,XX(LGHOND),MAXG,
     *              IA,PA,PB,FA,FB,FM2A,FM2B,XX(LDSH),NV)
C
         CALL RETFM(NEED)
      ELSE
         TDSKWRK = DSKWRK
         DSKWRK  = .TRUE.
         CALL SEQREW(IS)
         IF(LCFLAG) CALL SEQREW(LRFILE)
         CALL UAPMBDISK(SCFTYP,PA,PB,NBF,NBF2,FA,FB,FM2A,FM2B,
     *        BUF,IBUF,NINTMX,NOPK,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)
         IF(IPF.EQ.1) THEN
            CALL DDI_GSUMF(2313,FA,NV*NBF2)
            CALL DDI_GSUMF(2313,FB,NV*NBF2)
         END IF
         IF(IMF.EQ.1) THEN
            CALL DDI_GSUMF(2314,FM2A,NV*NBF3)
            CALL DDI_GSUMF(2314,FM2B,NV*NBF3)
         END IF
      END IF
C
      RETURN
      END
C
C*MODULE UTDDFT  *DECK UAPMBDISK
      SUBROUTINE UAPMBDISK(TYPSCF,PA,PB,L1,L2,FA,FB,FM2A,FM2B,
     *     XX,IX,NINTMX,NOPK,NV)
C
C     CSCAL : PARAMETER FOR COULOMB INTEGRALS
C     (CSCAL=1 FOR SINGLET AND CSCAL=0 FOR TRIPLET EXCITATIONS)
C     IMF : WHETHER TO CALCULATE (A-B) MATRIX OR NOT.
C     IPF : WHETHER TO CALCULATE (A+B) MATRIX OR NOT.
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION FA(L2,NV),FB(L2,NV),XX(NINTMX),IX(*)
      DIMENSION FM2A(L1,L1,NV),FM2B(L1,L1,NV)
      DIMENSION PA(L1,L1,NV),PB(L1,L1,NV)
C
      LOGICAL   LCFLAG,LRINT,UROHF
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      COMMON /ORDOPT/ NORDER(7)
      COMMON /PCKLAB/ LABSIZ
      COMMON /TDSCL/ CSCAL,IMF,IPF
C
      DATA HALF/0.5D+00/
      DATA TWO/2.0D+00/
      DATA UHF,ROHF/8HUHF     ,8HROHF    /
C
      UROHF=TYPSCF.EQ.UHF .OR. TYPSCF.EQ.ROHF
      CALL SEQREW(IS)
      IF(LCFLAG) CALL SEQREW(LRFILE)
      HFSCAL=DFTTYP(3)
      L3=L1*L1
      I = 0
      J = 0
      K = 0
      L = 0
      NXX = 0
      IF(LCFLAG) LRINT=.FALSE.
      IF(NOPK.NE.1) THEN
      WRITE(IW,*)'NOPK.NE.1'
      WRITE(IW,*)'TDDFT CALCULATION DOES NOT SUPPORT THIS INTEGRAL TYPE'
      CALL ABRT
      ENDIF
C
C     ----- INTEGRALS ARE NOT IN SUPERMATRIX FORM (NOPK=.TRUE.) -----
C
  620 CONTINUE
C     CALL PREAD(IS,XX,IX,NXX,NINTMX)
      CSCALT=CSCAL
      IF(LCFLAG)THEN
         IF(LRINT)THEN
            HFSCAL =1.0D+00
            CSCALT =0.0D+00
            CALL PREAD(LRFILE,XX,IX,NXX,NINTMX)
         ELSE
            HFSCAL =0.0D+00
            CSCALT =CSCAL
            CALL PREAD(IS,XX,IX,NXX,NINTMX)
         END IF
      ELSE
         CALL PREAD(IS,XX,IX,NXX,NINTMX)
      END IF
      IF(NXX .EQ. 0) GO TO 720
      NINT = IABS(NXX)
      IF(NINT .GT. NINTMX) CALL ABRT
      DO 700 M = 1,NINT
C
                       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
                       I = IPACK
                       J = JPACK
                       K = KPACK
                       L = LPACK
C
         VAL = XX(M)
         NIJ = IA(I)+J
         NKL = IA(K)+L
         IF(NORDER(7) .NE. 1) GO TO 630
C
C           USING SQUARE CANONICAL INTEGRAL FILE
C
         IF(NKL .GT. NIJ) GO TO 700
            IF(I .EQ. J) VAL=VAL*HALF
            IF(K .EQ. L) VAL=VAL*HALF
            IF(NIJ .EQ. NKL) VAL=VAL*HALF
C
  630 CONTINUE
      NIK = IA(I)+K
      NIL = IA(I)+L
      IF(J .LT. K) GO TO 640
      NJK = IA(J)+K
      NJL = IA(J)+L
      GO TO 680
C
  640 NJK = IA(K)+J
      IF(J .LT. L) GO TO 660
      NJL = IA(J)+L
      GO TO 680
C
  660 NJL = IA(L)+J
  680 CONTINUE
      VAL4 = (VAL+VAL)+(VAL+VAL)
      VAL4 = VAL4*CSCALT
      VAL  = VAL *HFSCAL
C
      DO 690 IMO=1,NV
      IF (IPF.EQ.0) GOTO 685
C     A+B
      IF(.NOT.UROHF) THEN
         DKL=PA(K,L,IMO)+PA(L,K,IMO)
         DKL=TWO*DKL
         DIJ=PA(I,J,IMO)+PA(J,I,IMO)
         DIJ=TWO*DIJ
         FA(NIJ,IMO) = FA(NIJ,IMO)+VAL4*DKL
         FA(NKL,IMO) = FA(NKL,IMO)+VAL4*DIJ
      ELSE
         DKL=PA(K,L,IMO)+PA(L,K,IMO)+PB(K,L,IMO)+PB(L,K,IMO)
         DIJ=PA(I,J,IMO)+PA(J,I,IMO)+PB(I,J,IMO)+PB(J,I,IMO)
         FA(NIJ,IMO) = FA(NIJ,IMO)+VAL4*DKL
         FA(NKL,IMO) = FA(NKL,IMO)+VAL4*DIJ
         FB(NIJ,IMO) = FB(NIJ,IMO)+VAL4*DKL
         FB(NKL,IMO) = FB(NKL,IMO)+VAL4*DIJ
      END IF
      DJL=PA(J,L,IMO)+PA(L,J,IMO)
      DJL=TWO*DJL
      DJK=PA(J,K,IMO)+PA(K,J,IMO)
      DJK=TWO*DJK
      DIL=PA(I,L,IMO)+PA(L,I,IMO)
      DIL=TWO*DIL
      DIK=PA(I,K,IMO)+PA(K,I,IMO)
      DIK=TWO*DIK
      FA(NIK,IMO) = FA(NIK,IMO)-VAL*DJL
      FA(NIL,IMO) = FA(NIL,IMO)-VAL*DJK
      FA(NJK,IMO) = FA(NJK,IMO)-VAL*DIL
      FA(NJL,IMO) = FA(NJL,IMO)-VAL*DIK
      IF(UROHF) THEN
         DJL=PB(J,L,IMO)+PB(L,J,IMO)
         DJL=TWO*DJL
         DJK=PB(J,K,IMO)+PB(K,J,IMO)
         DJK=TWO*DJK
         DIL=PB(I,L,IMO)+PB(L,I,IMO)
         DIL=TWO*DIL
         DIK=PB(I,K,IMO)+PB(K,I,IMO)
         DIK=TWO*DIK
         FB(NIK,IMO) = FB(NIK,IMO)-VAL*DJL
         FB(NIL,IMO) = FB(NIL,IMO)-VAL*DJK
         FB(NJK,IMO) = FB(NJK,IMO)-VAL*DIL
         FB(NJL,IMO) = FB(NJL,IMO)-VAL*DIK
      END IF
  685 CONTINUE
      IF (IMF.EQ.0) GOTO 690
C     A-B
      DJLIK=PA(L,J,IMO)-PA(J,L,IMO)
      DJLIK=TWO*DJLIK
      DJKIL=PA(K,J,IMO)-PA(J,K,IMO)
      DJKIL=TWO*DJKIL
      DILJK=PA(L,I,IMO)-PA(I,L,IMO)
      DILJK=TWO*DILJK
      DIKJL=PA(K,I,IMO)-PA(I,K,IMO)
      DIKJL=TWO*DIKJL
      FM2A(I,K,IMO) = FM2A(I,K,IMO)+VAL*DJLIK
      FM2A(K,I,IMO) = FM2A(K,I,IMO)-VAL*DJLIK
      FM2A(I,L,IMO) = FM2A(I,L,IMO)+VAL*DJKIL
      FM2A(L,I,IMO) = FM2A(L,I,IMO)-VAL*DJKIL
      FM2A(J,K,IMO) = FM2A(J,K,IMO)+VAL*DILJK
      FM2A(K,J,IMO) = FM2A(K,J,IMO)-VAL*DILJK
      FM2A(J,L,IMO) = FM2A(J,L,IMO)+VAL*DIKJL
      FM2A(L,J,IMO) = FM2A(L,J,IMO)-VAL*DIKJL
      IF(UROHF) THEN
         DJLIK=PB(L,J,IMO)-PB(J,L,IMO)
         DJLIK=TWO*DJLIK
         DJKIL=PB(K,J,IMO)-PB(J,K,IMO)
         DJKIL=TWO*DJKIL
         DILJK=PB(L,I,IMO)-PB(I,L,IMO)
         DILJK=TWO*DILJK
         DIKJL=PB(K,I,IMO)-PB(I,K,IMO)
         DIKJL=TWO*DIKJL
         FM2B(I,K,IMO) = FM2B(I,K,IMO)+VAL*DJLIK
         FM2B(K,I,IMO) = FM2B(K,I,IMO)-VAL*DJLIK
         FM2B(I,L,IMO) = FM2B(I,L,IMO)+VAL*DJKIL
         FM2B(L,I,IMO) = FM2B(L,I,IMO)-VAL*DJKIL
         FM2B(J,K,IMO) = FM2B(J,K,IMO)+VAL*DILJK
         FM2B(K,J,IMO) = FM2B(K,J,IMO)-VAL*DILJK
         FM2B(J,L,IMO) = FM2B(J,L,IMO)+VAL*DIKJL
         FM2B(L,J,IMO) = FM2B(L,J,IMO)-VAL*DIKJL
      END IF
  690 CONTINUE
  700 CONTINUE
C
      IF(NXX .GT. 0) GO TO 620
CC     DONE CONSTRUCTING MATRIX, OFF DIAGONAL ELEMENTS NEED HALVING
C
  720 CONTINUE
      IF(LCFLAG) THEN
         IF(LRINT) THEN
            LRINT=.FALSE.
         ELSE
            I = 0
            J = 0
            K = 0
            L = 0
            NXX  = 0
            LRINT=.TRUE.
            GO TO 620
         END IF
      END IF
C
      IF(IPF.EQ.1) THEN
         CALL DSCAL(L2*NV,HALF,FA,1)
         DO IMO=1,NV
            II = 0
            DO I=1,NUM
               II = II + I
               FA(II,IMO) = FA(II,IMO)+FA(II,IMO)
            END DO
         END DO
         IF (UROHF) THEN
            CALL DSCAL(L2*NV,HALF,FB,1)
            DO IMO=1,NV
               II = 0
               DO I=1,NUM
                  II = II + I
                  FB(II,IMO) = FB(II,IMO)+FB(II,IMO)
               END DO
            END DO
         END IF
      ENDIF
C
      IF(IMF.EQ.1) THEN
         CALL DSCAL(L3*NV,HALF,FM2A,1)
         DO IMO=1,NV
            DO I=1,NUM
               FM2A(I,I,IMO) = FM2A(I,I,IMO)*TWO
            END DO
         END DO
         IF (UROHF) THEN
            CALL DSCAL(L3*NV,HALF,FM2B,1)
            DO IMO=1,NV
               DO I=1,NUM
                  FM2B(I,I,IMO) = FM2B(I,I,IMO)*TWO
               END DO
            END DO
         END IF
      ENDIF
C
      CALL SEQREW(IS)
      IF(LCFLAG)CALL SEQREW(LRFILE)
      RETURN
C     ******
      END
C
C*MODULE UTDDFT  *DECK UTDTWOEI
      SUBROUTINE UTDTWOEI(TYPSCF,SCHWRZ,NINT,NSCHWZ,L1,L2,
     *                    XINTS,NSH2,GHONDO,MAXG,
     *                    IA,DA,DB,FA,FB,FM2A,FM2B,DSH,NV)
C
C     THIS IS A DUMMY CODE OF 'SUBROUTINE TWOEI' @ INT2A.SRC
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL CMBDIR
      LOGICAL OUT,SCHWRZ,SCHSKP,GOPARR,DSKWRK,MASWRK,DLB,SLB,C1GRP
      LOGICAL SKIPA,SKIPB,SKIPC,NPSYM
      LOGICAL PK,PANDK,BLOCK,GPSAVE
      LOGICAL LTRMST,LCFLAG,LRINT
      LOGICAL UROHF
C
      CHARACTER*8 INAOFL
C
      DIMENSION XINTS(NSH2),GHONDO(MAXG),IA(L1),DA(*),DB(*),
     *          FA(L2,*),FB(L2,*),FM2A(L1,L1,*),FM2B(L1,L1,*),DSH(NSH2)
      DIMENSION MI(48),MJ(48),MK(48),M0(48)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
      COMMON /ELGFIL/ NFILE,INAOFL
      COMMON /ELGTRM/ LTRMST,NFLTRM,NRCTRM,NPSTRM,NHTSHL
      COMMON /INT2IC/ NINTIC,ININTIC,NXXIC,LBUFPIC,LIXIC,LABSIX,NINTIX
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PKFIL / PK,PANDK,BLOCK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /SHLEXC/ NORGSH(3),NORGSP(3),IEXCH,NANGM,NGTH(4)
      COMMON /SHLG70/ ISH,JSH,KSH,LSH,IJKLXX(4)
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJ,IJ,KL,IJKL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TDSCL/ CSCAL,IMF,IPF
C
      PARAMETER (ZERO=0.0D+00)
      PARAMETER (HALF=0.5D+00)
      DATA UHF,ROHF/8HUHF     ,8HROHF    /
C
C          ----- TWO-ELECTRON INTEGRALS -----
C     THIS VERSION CAN HANDLE S,P,D,F,G AND L SHELLS
C
      UROHF=TYPSCF.EQ.UHF .OR. TYPSCF.EQ.ROHF
      TIM = ZERO
      CALL TSECND(TIM)
C
      CMBDIR= .TRUE.
C
C        THE OLD FASHIONED PARALLEL INTEGRAL TRANSFORMATIONS DO NOT
C        ALLOW THE AO INTEGRAL WORK TO BE RUN IN PARALLEL.  THE MODERN
C        DISTRIBUTED MEMORY TRANSFORMATIONS DO NOT PASS THROUGH HERE.
C
      GPSAVE = GOPARR
C
C     ----- INITIALIZATION FOR PARALLEL WORK -----
C     BOTH STATIC AND DYNAMIC LOAD BALANCING ARE IMPLEMENTED BELOW
C
      SLB = GOPARR  .AND.  IBTYP.EQ.0
      DLB = GOPARR  .AND.  IBTYP.EQ.1
      NEXT = -1
      MINE = -1
      IPCOUNT = ME - 1
C
      C1GRP = NT.EQ.1
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
      NGTH(4) = 1
      NGTH(3) = NGTH(4) * NANGM
      NGTH(2) = NGTH(3) * NANGM
      NGTH(1) = NGTH(2) * NANGM
      IF(NOPK.EQ.0) THEN
         NORGSH(1) = 0
         NORGSH(2) = NORGSH(1) + NANGM**4
         NORGSH(3) = NORGSH(2) + NANGM**4
         NORGSP(1) = 0
         NORGSP(2) = 256
         NORGSP(3) = 512
      ELSE
         DO I=1,3
            NORGSH(I) = 0
            NORGSP(I) = 0
         ENDDO
      END IF
C
      NINT  = 0
      NSCHWZ= 0
      SCHSKP=.FALSE.
      DENMAX = ZERO
C
C        NOW WE ARE READY TO LOOP OVER ALL NSHELL**4 SHELL QUARTETS
C
C     ----- I SHELL -----
C
      DO 920 II = IST,NSHELL
C
C     ----- CHECK CPU TIME -----
C
      CALL TSECND(TIM)
      IF(TIM.GE.TIMLIM) THEN
C        NOTHING CAN BE DONE FOR IN-CORE INTEGRALS: JUST FORGET THEM
         IF(MASWRK) WRITE(IW,9030)
         CALL ABRT
      END IF
C
C         ELONGATION METHOD
      IF(LTRMST) THEN
         IF(II.EQ.NHTSHL) THEN
            NFLTRM = NFILE
            NRCTRM = NREC
            NPSTRM = ICOUNT
         ENDIF
      ENDIF
C
C     ----- PRINT INTERMEDIATE RESTART DATA -----
C
      IF(NPRINT.NE.-5  .AND.  .NOT.CMBDIR .AND. MASWRK) THEN
         IF(ICOUNT.LE.NINTIC) THEN
            WRITE(IW,9015) II,JST,KST,LST,ICOUNT
         ELSE
            WRITE(IW,9010) II,JST,KST,LST,NREC,ICOUNT-NINTIC
         ENDIF
      ENDIF
C
C     ----- SKIP I SHELL IF NOT SYMMETRY UNIQUE -----
C     THIS, AND THE SIMILAR BRANCHINGS FOR THE J, K, AND L LOOPS IS
C     WHAT GENERATES THE "PETITE" RATHER THAN "GRANDE" INTEGRAL LIST.
C
      IF(C1GRP) THEN
         MI(1)=II
      ELSE
         DO 120 IT = 1,NT
            ID = MAPSHL(II,IT)
            IF (ID .GT. II) GO TO 920
            MI(IT) = ID
  120    CONTINUE
      END IF
C
C     ----- J SHELL -----
C
      J0 = JST
      DO 900 JJ = J0,II
      JST = 1
C
      IF(C1GRP) THEN
         MJ(1)=JJ
      ELSE
         DO 200 IT = 1,NT
            ID = MI(IT)
            JD = MAPSHL(JJ,IT)
            MJ(IT) = JD
            IF (ID .GE. JD) GO TO 160
            ND = ID
            ID = JD
            JD = ND
  160       IF (ID-II) 200,180,900
  180       IF (JD-JJ) 200,200,900
  200    CONTINUE
      END IF
C
C     ----- GO PARALLEL! -----
C
      IF (DLB) THEN
         MINE = MINE + 1
         IF (MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
         IF (NEXT.NE.MINE) GO TO 900
      END IF
C
C     ----- K SHELL -----
C
      K0 = KST
      DO 880 KK = K0,JJ
      KST = 1
C
      IF(C1GRP) THEN
         MK(1)=KK
      ELSE
         DO 340 IT = 1,NT
            ID = MI(IT)
            JD = MJ(IT)
            KD = MAPSHL(KK,IT)
            MK(IT) = KD
  240       IF (ID .GE. JD) GO TO 260
            ND = ID
            ID = JD
            JD = ND
  260       IF (JD .GE. KD) GO TO 280
            ND = JD
            JD = KD
            KD = ND
            GO TO 240
  280       IF (ID-II) 340,300,880
  300       IF (JD-JJ) 340,320,880
  320       IF (KD-KK) 340,340,880
  340    CONTINUE
      END IF
C
C     ----- GO PARALLEL! -----
C
      IF(SLB) THEN
         IPCOUNT = IPCOUNT + 1
         IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 880
      END IF
C
C     ----- L SHELL ----
C
      L0 = LST
      DO 860 LL = L0,KK
      LST = 1
C
      IF(C1GRP) THEN
         M0(1)=1
         N4=1
      ELSE
         N4 = 0
         DO 540 IT = 1,NT
            ID = MI(IT)
            JD = MJ(IT)
            KD = MK(IT)
            LD = MAPSHL(LL,IT)
  380       IF (ID .GE. JD) GO TO 400
            ND = ID
            ID = JD
            JD = ND
  400       IF (JD .GE. KD) GO TO 420
            ND = JD
            JD = KD
            KD = ND
            GO TO 380
  420       IF (KD .GE. LD) GO TO 440
            ND = KD
            KD = LD
            LD = ND
            GO TO 400
  440       IF (ID-II) 540,460,860
  460       IF (JD-JJ) 540,480,860
  480       IF (KD-KK) 540,500,860
  500       IF (LD-LL) 540,520,860
  520       N4 = N4+1
            M0(N4) = IT
  540    CONTINUE
      END IF
C
C         THE LOOP STRUCTURE IN THIS ROUTINE IS DESIGNED TO FACILITATE
C         SUPERMATRIX CONSTRUCTION BY HAVING UP TO THREE "EXCHANGED"
C         QUARTETS AVAILABLE AT ONCE.  THE LOOP STRUCTURE TO GENERATE
C         A MORE NORMAL CANONICAL ORDERING OF THE QUARTETS HITS THE
C         SAME QUARTETS IN A SLIGHTLY DIFFERENT ORDER, BUT BOTH LOOPS
C         WILL DO EXACTLY THE SAME QUARTETS.
C
C             CANONICAL                      SUPERMATRIX
C         DO ISH=1,NSHELL                 DO II=1,NSHELL
C           DO JSH=1,ISH                    DO JJ=1,II
C             IJSH = IA(ISH)+JSH
C             DO KSH=1,ISH                    DO KK=1,JJ
C               DO LSH=1,KSH                    DO LL=1,KK
C                 KLSH=IA(KSH)+LSH
C                 IF(IJSH.LT.KLSH),               [II JJ|KK LL],
C                    CYCLE KSH LOOP               [II KK|JJ LL],
C                 [ISH JSH|KSH LSH]               [II LL|JJ KK]
C               ENDDO                           ENDDO
C             ENDDO                           ENDDO
C           ENDDO                           ENDDO
C         ENDDO                           ENDDO
C
C     ----- CHECK FOR REDUNDANIES BETWEEN THE 3 COMBINATIONS -----
C            (II,JJ//KK,LL), (II,KK//JJ,LL), (II,LL//JJ,KK)
C
      SKIPA =  JJ.EQ.KK
      SKIPB = (II.EQ.KK) .OR. (JJ.EQ.LL)
      SKIPC = (II.EQ.JJ) .OR. (KK.EQ.LL)
      NPSYM = .FALSE.
      IF (SKIPA .OR. SKIPB .OR. SKIPC) GO TO 720
      NPSYM = .TRUE.
      DO 640 M = 1,N4
         IT = M0(M)
         IH = MI(IT)
         JH = MJ(IT)
         IF(JH.LE.IH) THEN
            ID = IH
            JD = JH
         ELSE
            ID = JH
            JD = IH
         END IF
         IF(.NOT.SKIPA) SKIPA = (ID.EQ.II .AND. JD.EQ.KK) .OR.
     *                          (ID.EQ.JJ .AND. JD.EQ.LL)
         IF(.NOT.SKIPB) SKIPB = (ID.EQ.II .AND. JD.EQ.LL) .OR.
     *                          (ID.EQ.JJ .AND. JD.EQ.KK)
         IF (SKIPA .AND. SKIPB) GO TO 660
         KH = MK(IT)
         IF(KH.LE.IH) THEN
            ID = IH
            KD = KH
         ELSE
            ID = KH
            KD = IH
         END IF
         IF(.NOT.SKIPC) SKIPC = (ID.EQ.II .AND. KD.EQ.LL) .OR.
     *                          (ID.EQ.JJ .AND. KD.EQ.KK)
         IF(SKIPA .AND. SKIPC) GO TO 680
         IF(SKIPB .AND. SKIPC) GO TO 700
  640 CONTINUE
      GO TO 720
C
  660 SKIPC = .TRUE.
      GO TO 720
  680 SKIPB = .TRUE.
      GO TO 720
  700 SKIPA = .TRUE.
C
C        GENERATE SYMMETRY FACTOR -Q4- FOR THIS QUARTET IN PETITE LIST
C
  720 CONTINUE
      Q4 = NT
      Q4 = Q4 / N4
C
C     ----- (II,JJ//KK,LL) -----
C
      IEXCH = 1
      ISH = II
      JSH = JJ
      KSH = KK
      LSH = LL
      QQ4 = Q4
      IF(SKIPA .AND. NPSYM) QQ4 = QQ4+Q4
      IF(SKIPB .AND. NPSYM) QQ4 = QQ4+Q4
      GO TO 780
C
C     ----- (II,KK//JJ,LL) -----
C
  740 IF (SKIPA) GO TO 760
      IEXCH = 2
      ISH = II
      JSH = KK
      KSH = JJ
      LSH = LL
      QQ4 = Q4
      IF (SKIPC .AND. NPSYM) QQ4 = QQ4+Q4
      GO TO 780
C
C     ----- (II,LL//JJ,KK) -----
C
  760 IF (SKIPB .OR. SKIPC) GO TO 840
      IEXCH = 3
      ISH = II
      JSH = LL
      KSH = JJ
      LSH = KK
      QQ4 = Q4
C
C        ----- COMPUTE TWO-ELECTRON INTEGRALS ----
C
  780 CONTINUE
      IF(PK .AND. IEXCH.EQ.1) CALL ZPKOUT(ISH,JSH,KSH,LSH,GHONDO,
     *                                    SKIPA,SKIPB,SKIPC,NPSYM)
C
C     APPLY THE SCHWARZ INEQUALITY TO SCREEN OUT SMALL INTEGRALS,
C      (II,JJ//KK,LL) .LE.  SQRT( (II,JJ//II,JJ)*(KK,LL//KK,LL) )
C     SEE, FOR EXAMPLE, J.L.WHITTEN, J.CHEM.PHYS. 58,4496-4501(1973)
C
      IF(SCHWRZ) THEN
         IJIJ = (ISH*ISH-ISH)/2 + JSH
         KLKL = (KSH*KSH-KSH)/2 + LSH
         TEST = QQ4*XINTS(IJIJ)*XINTS(KLKL)
         DENMAX = SCHWDN(DSH,ISH,JSH,KSH,LSH,IA)
         TEST = TEST*DENMAX
         SCHSKP = TEST.LT.CUTOFF
         IF(SCHSKP) NSCHWZ = NSCHWZ + 1
      END IF
      IF(SCHSKP) GO TO 820
C
C        ----- ELECTRON REPULSION INTEGRAL CALCULATION -----
C     THIS MAY USE ROTATED AXIS, ERIC, OR RYS QUADRATURE METHODS
C
      CALL SHELLQUART(ISH,JSH,KSH,LSH,GHONDO)
C
C        DIRECT-TDDFT
C
      CALL UTDDIRFCK(UROHF,DA,DB,FA,FB,FM2A,FM2B,
     *     GHONDO,L1,L2,NINT,NV,IA)
C
  820 CONTINUE
      GO TO (740,760,840),IEXCH
C
  840 CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
  860 CONTINUE
  880 CONTINUE
  900 CONTINUE
  920 CONTINUE
C
      IF(DLB) CALL DDI_DLBRESET
C
      IF(.NOT.LRINT) THEN
         IF(IPF.EQ.1) THEN
            CALL DSCAL(L2*NV,HALF,FA,1)
            DO IMO=1,NV
               II = 0
               DO I=1,L1
                  II = II + I
                  FA(II,IMO) = FA(II,IMO)+FA(II,IMO)
               END DO
            END DO
            IF (UROHF) THEN
               CALL DSCAL(L2*NV,HALF,FB,1)
               DO IMO=1,NV
                  II = 0
                  DO I=1,L1
                     II = II + I
                     FB(II,IMO) = FB(II,IMO)+FB(II,IMO)
                  END DO
               END DO
            END IF
         ENDIF
C
         IF(IMF.EQ.1) THEN
            CALL DSCAL(L1*L1*NV,HALF,FM2A,1)
            DO IMO=1,NV
               DO I=1,L1
                  FM2A(I,I,IMO) = FM2A(I,I,IMO)+FM2A(I,I,IMO)
               END DO
            END DO
            IF(UROHF) THEN
               CALL DSCAL(L1*L1*NV,HALF,FM2B,1)
               DO IMO=1,NV
                  DO I=1,L1
                     FM2B(I,I,IMO) = FM2B(I,I,IMO)+FM2B(I,I,IMO)
                  END DO
               END DO
            END IF
         ENDIF
C
      ENDIF
      GOPARR = GPSAVE
      RETURN
C
 9010 FORMAT(1X,'II,JST,KST,LST =',4I3,' NREC =',I10,' INTLOC =',I5)
 9015 FORMAT(1X,'II,JST,KST,LST =',4I3,' IN CORE, INTLOC =',I12)
 9030 FORMAT(//1X,'*** THIS JOB HAS EXHAUSTED ITS CPU TIME ***'/
     *         1X,'     (WHILE COMPUTING 2E- INTEGRALS)'///)
      END
C
C*MODULE UTDDFT  *DECK UTDDIRFCK
      SUBROUTINE UTDDIRFCK(UROHF,PA,PB,FA,FB,FM2A,FM2B,
     *                     GHONDO,NBF,NBF2,NINT,NV,IA)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL IANDJ,KANDL,SAME,OUT,LCFLAG,LRINT,UROHF
C
      DIMENSION GHONDO(*),IA(*)
      DIMENSION PA(NBF,NBF,NV),PB(NBF,NBF,NV)
      DIMENSION FA(NBF2,NV),FM2A(NBF,NBF,NV)
      DIMENSION FB(NBF2,NV),FM2B(NBF,NBF,NV)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /ERIOUT/ ISH,JSH,KSH,LSH,LSTRI,LSTRJ,LSTRK,LSTRL
      COMMON /MISC  / IANDJ,KANDL,SAME
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJX,IJ,KL,IJKL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
      COMMON /TDSCL/ CSCAL,IMF,IPF
C
      DATA HALF /0.5D+00/
      DATA TWO  /2.0D+00/
C
C-NEXT STATEMENT IS FOR VARIOUS IBM XLF 3.X AND 5.X COMPILERS-
C
      SAVE IJN,KLN
      HFSCAL=DFTTYP(3)
      CSCALT=CSCAL
      IF(LCFLAG) THEN
       IF(LRINT) THEN
       HFSCAL=1.0D+00
       CSCALT=0.0D+00
       ELSE
       HFSCAL=0.0D+00
       CSCALT=CSCAL
       ENDIF
      ENDIF
C     WRITE(6,*)'HFSCAL,CSCAL,IMF,NV=',HFSCAL,CSCAL,IMF,NV
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
C      WE NOW HAVE EVERYTHING READY, PUT INTEGRALS INTO FOCK MATRIX
C
                  NIJ = IA(II)+JJ
                  NKL = IA(KK)+LL
                  NIK = IA(II)+KK
                  NIL = IA(II)+LL
                  NJK = IA(JJ)+KK
                  NJL = IA(JJ)+LL
                  IF(JJ.LT.KK) NJK = IA(KK)+JJ
                  IF(JJ.LT.LL) NJL = IA(LL)+JJ
C
                  VAL2 = VAL+VAL
                  VAL4 = (VAL2+VAL2)*CSCALT
                  VAL  = VAL*HFSCAL
C
C     FOCK-LIKE MATRICES
C     A+B
                  DO 200 IMO=1,NV
                  IF (IPF.EQ.0) GOTO 190
C                 COULOMB
                  IF(.NOT.UROHF) THEN
                     DKL=PA(KK,LL,IMO)+PA(LL,KK,IMO)
                     DKL=TWO*DKL
                     DIJ=PA(II,JJ,IMO)+PA(JJ,II,IMO)
                     DIJ=TWO*DIJ
                     FA(NIJ,IMO) = FA(NIJ,IMO)+VAL4*DKL
                     FA(NKL,IMO) = FA(NKL,IMO)+VAL4*DIJ
                  ELSE
                     DKL=PA(KK,LL,IMO)+PA(LL,KK,IMO)+
     *                   PB(KK,LL,IMO)+PB(LL,KK,IMO)
                     DIJ=PA(II,JJ,IMO)+PA(JJ,II,IMO)+
     *                   PB(II,JJ,IMO)+PB(JJ,II,IMO)
                     FA(NIJ,IMO) = FA(NIJ,IMO)+VAL4*DKL
                     FA(NKL,IMO) = FA(NKL,IMO)+VAL4*DIJ
                     FB(NIJ,IMO) = FB(NIJ,IMO)+VAL4*DKL
                     FB(NKL,IMO) = FB(NKL,IMO)+VAL4*DIJ
                  END IF
C                 EXCHANGE
                  DJL=PA(JJ,LL,IMO)+PA(LL,JJ,IMO)
                  DJL=TWO*DJL
                  DJK=PA(JJ,KK,IMO)+PA(KK,JJ,IMO)
                  DJK=TWO*DJK
                  DIL=PA(II,LL,IMO)+PA(LL,II,IMO)
                  DIL=TWO*DIL
                  DIK=PA(II,KK,IMO)+PA(KK,II,IMO)
                  DIK=TWO*DIK
                  FA(NIK,IMO) = FA(NIK,IMO)-VAL*DJL
                  FA(NIL,IMO) = FA(NIL,IMO)-VAL*DJK
                  FA(NJK,IMO) = FA(NJK,IMO)-VAL*DIL
                  FA(NJL,IMO) = FA(NJL,IMO)-VAL*DIK
                  IF (UROHF) THEN
                     DJL=PB(JJ,LL,IMO)+PB(LL,JJ,IMO)
                     DJL=TWO*DJL
                     DJK=PB(JJ,KK,IMO)+PB(KK,JJ,IMO)
                     DJK=TWO*DJK
                     DIL=PB(II,LL,IMO)+PB(LL,II,IMO)
                     DIL=TWO*DIL
                     DIK=PB(II,KK,IMO)+PB(KK,II,IMO)
                     DIK=TWO*DIK
                     FB(NIK,IMO) = FB(NIK,IMO)-VAL*DJL
                     FB(NIL,IMO) = FB(NIL,IMO)-VAL*DJK
                     FB(NJK,IMO) = FB(NJK,IMO)-VAL*DIL
                     FB(NJL,IMO) = FB(NJL,IMO)-VAL*DIK
                  END IF
  190             CONTINUE
                  IF (IMF.EQ.0) GOTO 200
C     A-B
                  DJLIK=PA(LL,JJ,IMO)-PA(JJ,LL,IMO)
                  DJLIK=TWO*DJLIK
                  DJKIL=PA(KK,JJ,IMO)-PA(JJ,KK,IMO)
                  DJKIL=TWO*DJKIL
                  DILJK=PA(LL,II,IMO)-PA(II,LL,IMO)
                  DILJK=TWO*DILJK
                  DIKJL=PA(KK,II,IMO)-PA(II,KK,IMO)
                  DIKJL=TWO*DIKJL
                  FM2A(II,KK,IMO) = FM2A(II,KK,IMO)+VAL*DJLIK
                  FM2A(KK,II,IMO) = FM2A(KK,II,IMO)-VAL*DJLIK
                  FM2A(II,LL,IMO) = FM2A(II,LL,IMO)+VAL*DJKIL
                  FM2A(LL,II,IMO) = FM2A(LL,II,IMO)-VAL*DJKIL
                  FM2A(JJ,KK,IMO) = FM2A(JJ,KK,IMO)+VAL*DILJK
                  FM2A(KK,JJ,IMO) = FM2A(KK,JJ,IMO)-VAL*DILJK
                  FM2A(JJ,LL,IMO) = FM2A(JJ,LL,IMO)+VAL*DIKJL
                  FM2A(LL,JJ,IMO) = FM2A(LL,JJ,IMO)-VAL*DIKJL
C
                  IF (UROHF) THEN
                     DJLIK=PB(LL,JJ,IMO)-PB(JJ,LL,IMO)
                     DJLIK=TWO*DJLIK
                     DJKIL=PB(KK,JJ,IMO)-PB(JJ,KK,IMO)
                     DJKIL=TWO*DJKIL
                     DILJK=PB(LL,II,IMO)-PB(II,LL,IMO)
                     DILJK=TWO*DILJK
                     DIKJL=PB(KK,II,IMO)-PB(II,KK,IMO)
                     DIKJL=TWO*DIKJL
                     FM2B(II,KK,IMO) = FM2B(II,KK,IMO)+VAL*DJLIK
                     FM2B(KK,II,IMO) = FM2B(KK,II,IMO)-VAL*DJLIK
                     FM2B(II,LL,IMO) = FM2B(II,LL,IMO)+VAL*DJKIL
                     FM2B(LL,II,IMO) = FM2B(LL,II,IMO)-VAL*DJKIL
                     FM2B(JJ,KK,IMO) = FM2B(JJ,KK,IMO)+VAL*DILJK
                     FM2B(KK,JJ,IMO) = FM2B(KK,JJ,IMO)-VAL*DILJK
                     FM2B(JJ,LL,IMO) = FM2B(JJ,LL,IMO)+VAL*DIKJL
                     FM2B(LL,JJ,IMO) = FM2B(LL,JJ,IMO)-VAL*DIKJL
                  END IF
  200             CONTINUE
C
  300          CONTINUE
  320       CONTINUE
  340    CONTINUE
  360 CONTINUE
C
      RETURN
      END
