C 11 APR 08 - HPTI - GRIDS AND PARALLELIZATION
C 11 APR 08 - MC  - FIX SPHERICAL HARMONIC GRADIENT
C  4 MAR 08 - MWS - ADD EXETYP=CHECK MECHANISM
C 20 AUG 07 - MC  - NEW MODULE FOR TDDFT ENERGY GRADIENT COMPUTATION
C
C*MODULE TDDGRD  *DECK TDGRAD
      SUBROUTINE TDGRAD
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,GOPSAV,LCFLAG,LRINT,TRIPLET,SG1,SG1T
C
      PARAMETER (MXATM=2000, MXRT=100)
      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 /DFLEB/  NLEB(MXGRID),NLEB0(MXGRID)
      COMMON /FMCOM / X(1)
      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 /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
C     ----- TD-DFT RESPONSE EQUATION DRIVER -----
C     WRITTEN BY MAHITO CHIBA, NATIONAL INSTITUTE OF ADVANCED
C     INDUSTRIAL SCIENCE AND TECHNOLOGY (AIST), TSUKUBA, JAPAN.
C
C     PARALLEL VERSION OF THE GRADIENT CALCULATION ONLY WORKS FOR
C     DYANAMIC LOAD BALANCING
C
      GOPSAV=GOPARR
      NPROCS=NPROC
C
      IF(GOPARR) THEN
         IF(IBTYP.EQ.0) THEN
            WRITE(*,*)'STATIC LOAD BALANCING IS NOT SUPPORTED'
            WRITE(*,*)'IN THE PARALLEL TDDFT GRADIENT CALCULATION'
            WRITE(*,*)'DISABLING PARALLELIZATION'
            GOPARR=.FALSE.
            NPROC=1
            IF(.NOT.MASWRK) GOTO 100
         ENDIF
      ENDIF
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1 * L1
      IF(LCFLAG) THEN
         IF(MASWRK) WRITE(IW,1000)
      ELSE
         IF(MASWRK) WRITE(IW,1100)
      ENDIF
C
      IF(TRIPLET) THEN
        WRITE(IW,*) 'TRIPLET IS NOT SUPPORTED FOR TDDFT-GRADIENT'
        CALL ABRT
      END IF
C
      IF (NTHST.GT.MXRT) THEN
        WRITE(IW,*) 'TDDFT-GRADIENT CALCULATION IS ALLOWED ONLY',
     *              'FOR IROOT.LE.',MXRT
        CALL ABRT
      ENDIF
C
C     -- SWITCH TO THE SMALL GRID SIZE FOR TDDFT
C
      NRADS=NRAD
      NTHES=NTHE
      NPHIS=NPHI
      NLEBS=NLEB(1)
      NRAD=NRADT
      NTHE=NTHET
      NPHI=NPHIT
      NLEB(1)=NLEBT
C
C        POINT GROUP SYMMETRY CANNOT BE USED DURING THE RESPONSES.
C        THE SELECTED STATE MAY NOT HAVE A SYMMETRIC DENSITY, SO
C        WE ALSO HAVE TO DO THE GRADIENT INTEGRAL TERMS (LATER)
C        W/O SYMMETRY.
C
      CALL SYMOFF
      CALL RPADEN
      CALL SYMON
C
C     AFTER THIS, DAF #16 HOLDS THE FULL TDDFT DENSITY MATRIX,
C     WHEREAS #308 WILL CONTAIN THE GROUND STATE DENSITY.
C
      CALL VALFM(LOADFM)
      LDD  = LOADFM + 1
      LD   = LDD    + L3
      LD0  = LD     + L2
      LAST = LD0    + L2
      NEED = LAST - LOADFM -1
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) THEN
         CALL VCLR(X(LD) ,1,L2)
         CALL VCLR(X(LD0),1,L2)
      ELSE
         CALL DENDD1(X(LD),X(LD0),L2)
         CALL DCOPY(L2,X(LD),1,X(LD0),1)
         CALL DENTD(X(LD),X(LDD),L1)
      END IF
      CALL DAWRIT(IDAF,IODA,X(LD0),L2,308,0)
      CALL DAWRIT(IDAF,IODA,X(LD),L2,16,0)
      CALL RETFM(NEED)
C
      IF(MASWRK) WRITE(IW,9000)
      CALL TIMIT(1)
C
      IF(NPRINT.NE.-5) CALL PROPTY('TDDF')
C
C     -- SWITCH BACK TO THE REGULAR GRID SIZE, E.G. FOR GRADIENT STEP
C
      NRAD=NRADS
      NTHE=NTHES
      NPHI=NPHIS
      NLEB(1)=NLEBS
C
  100 CONTINUE
      NPROC =NPROCS
      GOPARR=GOPSAV
      RETURN
C
 1000 FORMAT(/1X,71(1H-)/
     *  19X,'LC-TDDFT ENERGY GRADIENT CALCULATION'//
     *  14X,'BY MAHITO CHIBA, TAKAO TSUNEDA, KIMIHIKO HIRAO,'/
     *  21X,'J. CHEM. PHYS. 124 (2006) 144106.'/
     *       1X,71(1H-)/)
 1100 FORMAT(/1X,71(1H-)/
     *  19X,'TD-DFT ENERGY GRADIENT CALCULATION'/
     *  19X,'CODED BY MAHITO CHIBA (AIST, JAPAN)'/
     *       1X,71(1H-)/)
 9000 FORMAT(1X,'..... DONE WITH TD-DFT EXCITED STATE RESPONSE',
     *          ' AND DENSITY MATRIX .....')
      END
C
C*MODULE TDDGRD  *DECK RPADEN
      SUBROUTINE RPADEN
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000,MXRT=100,ZERO=0.0D+00,ONE=1.0D+00)
      PARAMETER (MXGRID=10)
C
      LOGICAL TRIPLET,PACK2E,GOPARR,DSKWRK,MASWRK,SG1,SG1T
C
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /FMCOM / X(1)
      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 /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
      L0=NINTMX
      L1= NUM
      L2=(NUM*(NUM+1))/2
      L3= NUM* NUM
      LX=NQMT
      NOCC=NA
      NVIR=LX-NOCC
      L7=NOCC*NVIR
C
C     --- MEMORY ALLOCATION
C
      CALL VALFM(LOADFM)
      I10=LOADFM+1
C     I10:TA(L3)OR PA(L3)
      I20=I10+L3
C     I20:V1(L7)
      I30=I20+L7
C     I30:V2(L7)
      I40=I30+L7
C     I40:GXC(L2)
      I50=I40+L2
C     I50:FXCV(L2)
      I60=I50+L2
C     I60:FXCT(L2)
      I70=I60+L2
C     I70:R(L7) OR Z(L7)
      I80=I70+L7
C     I80:QAI(L7)
      I90=I80+L7
C     I90:TEMP
      I100=I90+L3
C     I100:MO(L3)
      I110=I100+L3
C     I110:E(LX)
      LAST=I110+LX
      NEEDRPA=LAST-I10
      CALL GETFM(NEEDRPA)
C
C     ---- READ MO
C
      CALL DAREAD(IDAF,IODA,X(I100),L1*LX,15,0)
C
C     ---- CALCULATE TA(L1,L1)
C
      IF(EXETYP.EQ.CHECK) GO TO 50
C
      CALL DAREAD(IDAF,IODA,X(I20),L7,IRECTD,0)
      CALL DAREAD(IDAF,IODA,X(I30),L7,IRECTD+1,0)
      CALL TCONST(X(I10),X(I20),X(I30),LX,NOCC)
C     SUBROUTINE TCONST(T,V1,V2,LX,NOCC)
      CALL DGEMM('N','N',L1,LX,LX,ONE,X(I100),L1,X(I10),LX,ZERO
     * ,X(I90),L1)
      CALL DGEMM('N','T',L1,L1,LX,ONE,X(I90),L1,X(I100),L1,ZERO
     * ,X(I10),L1)
C
C     ---- CALCULATE GXC(L2),FXCV(L2),FXCT(L2):GRID CALCULATION
C
   50 CONTINUE
      IF(NDFTFG.EQ.1) THEN
         CALL FGMNG(X(I40),X(I50),X(I60),X(I10),
     *              L1,L2,L3,LX,IRECTD,NEEDMNG)
      ENDIF
      NEEDMX = NEEDMNG + NEEDRPA
C
C     ---- R(L7)
C
      CALL RMNG(X(I70),X(I80),X(I40),X(I50),X(I60),X(I10),X(I20),
     *          X(I30),X(I100),NOCC,NVIR,L0,L1,L2,L3,L7,LX,
     *          LABSIZ,IRECTD,NEEDMNG)
      NEEDMX = MAX(NEEDMX,NEEDMNG + NEEDRPA)
C     SUBROUTINE RMNG(R,QAI,GXC,FXCV,FXCT,TA,V1,V2,V,NOCC,NVIR,L0,
C    * L1,L2,L3,L7,LX,LABSIZ,IRECTD)
C
C     ---- Z_MO(L7) AND PA(L1,L1) =TA+ZA
C
      CALL DAREAD(IDAF,IODA,X(I110),LX,17,0)
      CALL ZPMNG(X(I70),X(I10),X(I100),X(I110),L0,L1,L2,L3,LX,L7,
     *           NOCC,LABSIZ,IRECTD,NEEDMNG)
      NEEDMX = MAX(NEEDMX,NEEDMNG + NEEDRPA)
C
C     ---- CALCULATE W_AO(L1,L1)
C
      CALL WMNG(X(I80),X(I40),X(I10),X(I70),X(I20),X(I30),X(I100),
     *          X(I110),NOCC,L0,L1,L2,L3,LX,LABSIZ,NEEDMNG)
      NEEDMX = MAX(NEEDMX,NEEDMNG + NEEDRPA)
      IF(MASWRK) WRITE(IW,9010) NEEDMX
C
      CALL RETFM(NEEDRPA)
      RETURN
 9010 FORMAT(1X,'MAXIMUM MEMORY DURING TDDFT RESPONSES=',I10,' WORDS.')
      END
C*MODULE TDDGRD  *DECK TCONST
      SUBROUTINE TCONST(T,V1,V2,LX,NOCC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION T(LX,LX),V1(*),V2(*)
      DATA ZERO/0.0D+00/
      DATA ONE /1.0D+00/
      DATA TWO /2.0D+00/
      CALL VCLR(T,1,LX*LX)
C
C     -- T(VIR/VIR)
C
      FAC=ONE/TWO
      DO I=NOCC+1,LX
      DO J=NOCC+1,LX
      II=I-NOCC
      JJ=J-NOCC
        DUMK=ZERO
      DO K=1,NOCC
      KI=(II-1)*NOCC+K
      KJ=(JJ-1)*NOCC+K
      DUM=V1(KI)*V1(KJ)+V2(KI)*V2(KJ)
      DUMK=DUMK+DUM
      ENDDO
      T(I,J)=FAC*DUMK
      ENDDO
      ENDDO
C
C     -- T(OCC/OCC)
C
      FAC=-ONE/TWO
      DO I=1,NOCC
      DO J=1,NOCC
        DUMK=ZERO
      DO K=NOCC+1,LX
      KK=K-NOCC
      KI=(KK-1)*NOCC+I
      KJ=(KK-1)*NOCC+J
      DUM=V1(KI)*V1(KJ)+V2(KI)*V2(KJ)
      DUMK=DUMK+DUM
      ENDDO
      T(I,J)=FAC*DUMK
      ENDDO
      ENDDO
C
C     -- T(OCC/VIR) ELEMENTS ARE ZERO
C
      RETURN
      END
C*MODULE TDDGRD  *DECK FGMNG
      SUBROUTINE FGMNG(GXC,FXCV,FXCT,TA,L1,L2,L3,LX,IRECTD,NEED)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION GXC(*),FXCV(*),FXCT(*),TA(*)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /FMCOM / X(1)
      COMMON /INFGRD/RHOMIN,ILENG,MAXGRD
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
C     --- MEMORY ALLOCATION
C
      CALL VALFM(LOADFM)
      I10=LOADFM+1
C     I10:V1A(L3) OR TA(L3)
      LASTI=I10+L3
C     -- MEMORY FOR DFT GRID
      J10=LASTI+1
C     J10:XYZGRD(MAXGRD*3)
      J20=J10+MAXGRD*3
C     J20:XYZWGT(MAXGRD)
      J30=J20+MAXGRD
C     J30:XYZDCH
      J40=J30+4*L1*ILENG
C     J40:V(L3)ORD(L2)
      J50=J40+L3
C     J50:---
      J90=J50
C     J90:KCP(MAXGRD)
      J100=J90 +MAXGRD
C     J100:COEF(3*ILENG)
      J110=J100+3*ILENG
C     J110:DTRAI_(4*ILENG)
      J120=J110+4*ILENG
C     J120:DRHOI_(4*MAXGRD)
      J130=J120+4*MAXGRD
C     J130:EX0(ILENG)
      J140=J130+ILENG
C     J140:EX(18*ILENG)
      J150=J140+18*ILENG
C     J150:EC0(ILENG)
      J160=J150+ILENG
C     J160:EC(55*ILENG)
      J170=J160+55*ILENG
      J180=J170
C     J180:AOMAX(L1)
      J190=J180+L1
C     J190:IAO(L1)
      J200=J190+L1
C     J200:IPRGRD(2*L3)
      LAST=J200+2*L3
      IDUM=LAST
      NEED=LAST-I10
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C     ----- CONSTRUCT DENSITY
C     -- READ MO
      CALL DAREAD(IDAF,IODA,X(J40),L1*LX,15,0)
      CALL TDDFTSET(RHOSET,X(J10),X(J20),X(J30),X(J40),
     *              X(J120),X(J180),X(J190),ILENG,NPTGRD,L1)
C        ROUTINE TDDFTSET(RHOSET,XYZGRD,XYZWGT,XYZDCH,V,
C    *                    DRHOI_,AOMAX,GMO,ILENG,NPTGRD,L1)
C
C     ---- CALCULATE GXC(L2)
C                               --READ V1A
      CALL DAREAD(IDAF,IODA,X(I10),L3,IRECTD+2,0)
      CALL VCLR(GXC,1,L2)
      CALL TDFXCP(GXC,RHO,X(J10),X(J20),X(J30),X(I10),X(IDUM),
     * X(J120),X(J110),X(J100),X(J140),X(J160),X(J130),X(J150),X(J180),
     * X(J200),X(J190),ILENG,NPTGRD,L1,L2,LX,3,2)
C
C        ROUTINE TDFXCP(F,RHO,XYZGRD,XYZWGT,XYZDCH,VA,VDUM,
C    * DRHOI_,DTRAI_,COEF,EX,EC,EX0,EC0,AOMAX,VALGA,IAO,
C    * ILENG,NPTGRD,L1,L2,LX,NDER,IMETHOD)
C
      IF(GOPARR.AND.IBTYP.EQ.1) THEN
         CALL DDI_GSUMF(2332,GXC, L2)
         CALL DDI_GSUMF(2333,RHO,1)
      ENDIF
C
C     ---- CALCULATE FXCV(L2)
C
      CALL VCLR(FXCV,1,L2)
      CALL TDFXCP(FXCV,RHO,X(J10),X(J20),X(J30),X(I10),X(IDUM),
     * X(J120),X(J110),X(J100),X(J140),X(J160),X(J130),X(J150),X(J180),
     * X(J200),X(J190),ILENG,NPTGRD,L1,L2,LX,2,2)
      IF(GOPARR.AND.IBTYP.EQ.1) THEN
         CALL DDI_GSUMF(2334,FXCV, L2)
         CALL DDI_GSUMF(2335,RHO,1)
      ENDIF
C
C     ---- CALCULATE FXCT(L2)
C
      CALL VCLR(FXCT,1,L2)
      CALL TDFXCP(FXCT,RHO,X(J10),X(J20),X(J30),TA,X(IDUM),
     * X(J120),X(J110),X(J100),X(J140),X(J160),X(J130),X(J150),X(J180),
     * X(J200),X(J190),ILENG,NPTGRD,L1,L2,LX,2,2)
      IF(GOPARR.AND.IBTYP.EQ.1) THEN
         CALL DDI_GSUMF(2336,FXCT, L2)
         CALL DDI_GSUMF(2337,RHO,1)
      ENDIF
C
  800 CONTINUE
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE TDDGRD  *DECK RMNG
      SUBROUTINE RMNG(R,QAI,GXC,FXCV,FXCT,TA,V1,V2,V,NOCC,NVIR,L0,
     *                L1,L2,L3,L7,LX,LABSIZ,IRECTD,NEED)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION R(*),QAI(*),GXC(*),FXCV(*),FXCT(*),TA(*),
     *          V1(*),V2(*),V(*)
C
      LOGICAL SG1
C
      PARAMETER (MXGRID=10)
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
C
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
      CALL VALFM(LOADFM)
      CALL GOTFM(NGOTMX)
      I10 = LOADFM + 1
C     I10:V1A(L1,L1) OR SCRATCH
      I20=I10+L3
C     I20:SCRATCH
      I30=I20+L3
C     I30:---
      I40=I30
C     I40:XX(L0)
      I50=I40+L0
C     I50:IX(L0*LABSIZ)
      I60=I50+L0*LABSIZ
C     I60:FP1(L2)
      I70=I60+L2
C     I70:FPT(L2)
      I80=I70+L2
C     I80:FM2(L1,L1)
      I90=I80+L3
C     I90:--
      I100=I90
C     I100:--
      I110=I100
C     I110:--
      I120=I110
C     I120:HPV1IJ(NOCC,NOCC)
      I130=I120+NOCC*NOCC
C     I130:HPV1AB(NVIR,NVIR)
      I140=I130+NVIR*NVIR
C     I140:HMV2IJ(NOCC,NOCC)
      I150=I140+NOCC*NOCC
C     I150:HMV2AB(NVIR,NVIR)
      I155=I150+NVIR*NVIR
C     I155:HPTIA(L7)
      I160=I155+L7
C     I160:--
      I170=I160
C     I170:G_IA(L7)
      LAST=I170+L7
      IDUM=LAST
      NEED=LAST-I10
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C     ---- FPT(L2),FP1(L2),FM2(L1,L1)
C     FPT --TA
      CALL VCLR(X(I70),1,L2)
      CALL TD2E(TA,X(I70),X(IDUM),X(I40),X(I50),L1,ONE,0,1,1)
C     SUBROUTINE TD2E(PA,F,FM2,BUF,IBUF,NBF,SCAL,IMF,IPF,NV)
      IF(NDFTFG.EQ.1) CALL DAXPY(L2,2.0D+00,FXCT,1,X(I70),1)
C     FP1 --V1A
      CALL DAREAD(IDAF,IODA,X(I10),L3,IRECTD+2,0)
      CALL VCLR(X(I60),1,L2)
      CALL TD2E(X(I10),X(I60),X(IDUM),X(I40),X(I50),L1,ONE,0,1,1)
      IF(NDFTFG.EQ.1) CALL DAXPY(L2,2.0D+00,FXCV,1,X(I60),1)
C     FM2 --V2A
      CALL DAREAD(IDAF,IODA,X(I10),L3,IRECTD+3,0)
      CALL VCLR(X(I80),1,L3)
      CALL TD2E(X(I10),X(IDUM),X(I80),X(I40),X(I50),L1,ONE,1,0,1)
C
C     ---- HPV1AB, HMV2AB, HPV1IJ, HMV2IJ, HPTIA
C     --HPV1AB
      CALL CANTOGEN(X(I60),X(I10),L1)
      CALL MNTOAB(X(I10),X(I130),X(I20),V,L1,LX,NOCC,NVIR)
C     SUBROUTINE MNTOAB(PAO,PMO,SCR,V,L1,LX,NOCC,NVAL)
C     --HMV2AB
      CALL MNTOAB(X(I80),X(I150),X(I20),V,L1,LX,NOCC,NVIR)
C     --HPV1IJ
      CALL CANTOGEN(X(I60),X(I10),L1)
      CALL DGEMM('N','N',L1,NOCC,L1,ONE,X(I10),L1,V,L1,ZERO
     * ,X(I20),L1)
      CALL DGEMM('T','N',NOCC,NOCC,L1,ONE,V,L1,X(I20),L1,ZERO
     * ,X(I120),NOCC)
C     --HMV2IJ
      CALL DGEMM('N','N',L1,NOCC,L1,ONE,X(I80),L1,V,L1,ZERO
     * ,X(I10),L1)
      CALL DGEMM('T','N',NOCC,NOCC,L1,ONE,V,L1,X(I10),L1,ZERO
     * ,X(I140),NOCC)
C     --HPTIA
      CALL CANTOGEN(X(I70),X(I10),L1)
      CALL MNTOIA(X(I10),X(I155),V,X(I20),L1,LX,L7,NOCC,1)
C
C     ---- GIA FOR R (DFT)
C
      IF(NDFTFG.EQ.1) THEN
       CALL CANTOGEN(GXC,X(I10),L1)
       CALL MNTOIA(X(I10),X(I170),V,X(I20),L1,LX,L7,NOCC,1)
      ELSE
       CALL VCLR(X(I170),1,L7)
      ENDIF
C
C     ---- R(L7)
C
      CALL RCALC(R,QAI,V1,V2,X(I170),X(I120),
     * X(I140),X(I130),X(I150),X(I155),LX,L7,
     * NOCC,NOCC,NVIR)
C
  800 CONTINUE
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE TDDGRD  *DECK RCALC
      SUBROUTINE RCALC(R,QAI,V1,V2,GIA,HPV1IJ,HMV2IJ,HPV1AB,
     *                 HMV2AB,HPTIA,LX,L7,NOCC,NOCCG,NVALG)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION R(*),QAI(*),V1(*),V2(*),GIA(*),HPV1AB(NVALG,*),
     *          HMV2AB(NVALG,*),HPV1IJ(NOCCG,*),HMV2IJ(NOCCG,*),
     *          HPTIA(*)
C
      LOGICAL SG1
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
C
      DATA ZERO/0.0D+00/
C
      DO I=1,NOCC
      DO J=NOCC+1,LX
      JJ=J-NOCC
      IJ=(JJ-1)*NOCC+I
C     --- FOR DUM1
        DUM1=ZERO
        DO K=NOCC+1,LX
      KK=K-NOCC
      IB=(KK-1)*NOCC+I
        DUM=V1(IB)*HPV1AB(JJ,KK)+V2(IB)*HMV2AB(JJ,KK)
      DUM1=DUM1+DUM
      END DO
C     --- FOR DUM2
        DUM2=ZERO
        DO K=1,NOCC
      KA=(JJ-1)*NOCC+K
        DUM2=DUM2+V1(KA)*HPV1IJ(K,I)+V2(KA)*HMV2IJ(K,I)
      END DO
      QAI(IJ)=DUM2
C     --- FOR DUM3
      DUM3=HPTIA(IJ)
C     --- TOTAL
C       WRITE(6,*)'DUM1,DUM2,DUM3=',DUM1,DUM2,DUM3
      R(IJ)= DUM1-DUM2+DUM3
      END DO
      END DO
C     --- ADD GXC TO R
      IF(NDFTFG.EQ.1) CALL DAXPY(L7,2.0D+00,GIA,1,R,1)
C
      CALL DSCAL(L7,-1.0D+00,R,1)
C
      RETURN
      END
C*MODULE TDDGRD  *DECK MNTOAB
      SUBROUTINE MNTOAB(PAO,PMO,SCR,V,L1,LX,NOCC,NVAL)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (ZERO=0.0D+00)
      DIMENSION PAO(L1,*)
      DIMENSION PMO(NVAL,*)
      DIMENSION SCR(L1,*)
      DIMENSION V(L1,LX)
C
C     -- (M,N) TO (M,B)
C
      DO M=1,L1
      DO J=NOCC+1,LX
      JJ=J-NOCC
       DUMN=ZERO
       DO N=1,L1
       DUMN=DUMN+PAO(M,N)*V(N,J)
       END DO
      SCR(M,JJ)=DUMN
      END DO
      END DO
C
C     -- (M,B) TO (A,B)
C
      DO I=NOCC+1,LX
      DO J=NOCC+1,LX
      II=I-NOCC
      JJ=J-NOCC
       DUMM=ZERO
       DO M=1,L1
       DUMM=DUMM+SCR(M,JJ)*V(M,I)
       END DO
      PMO(II,JJ)=DUMM
      END DO
      END DO
      RETURN
      END
C*MODULE TDDGRD  *DECK ZPMNG
      SUBROUTINE ZPMNG(R,TA,V,E,L0,L1,L2,L3,LX,L7,NOCC,LABSIZ,
     *                 IRECTD,NEED)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SG1
C
      DIMENSION R(*),TA(*),V(*),E(*)
C
      PARAMETER (MXGRID=10)
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
C
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /FMCOM / X(1)
      COMMON /INFGRD/RHOMIN,ILENG,MAXGRD
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
      DATA THRE/1.0D-10/
      DATA MXITER/50/
C
      IF(MASWRK) THEN
      WRITE(IW,*)'   -------------------------'
      WRITE(IW,*)'      START Z-VECTOR LOOP'
      WRITE(IW,*)'   -------------------------'
      ENDIF
C
      CALL VALFM(LOADFM)
      CALL GOTFM(NGOTMX)
      I10=LOADFM+1
C     I10:--
      I20=I10
C     I20:XM(L7)
      I30=I20+L7
C     I30:XM_IN(L7)
      I40=I30+L7
C     I40:X_T(L7)
      K10=I40+L7
C     K10:PK_AO(L1,L1) OR SCRATCH
      K20=K10+L3
C     K20:AP_AO(L2)
      K30=K20+L2
C     K30:APXC(L2)
      K40=K30+L2
C     K40:XX(L0)
      K50=K40+L0
C     K50:IX(L0*LABSIZ)
      K60=K50+L0*LABSIZ
C     K60:---
      I50=K60
C     I50:A_PK(L7)
      I60=I50+L7
C     I60:---
      I70=I60
C     I70:R(L7)
      I80=I70+L7
C     I80:PK(L7)
      I90=I80+L7
      I100=I90
C     I100:Z_AO(L1,L1) OR SCRACTH
      LASTI=I100+L1*L1
      IF (NDFTFG.EQ.1) THEN
C
C     -- MEMORY FOR DFT GRID
C
      J10=LASTI+1
C     J10:XYZGRD(MAXGRD*3)
      J20=J10+MAXGRD*3
C     J20:XYZWGT(MAXGRD)
      J30=J20+MAXGRD
C     J30:XYZDCH
      J40=J30+4*L1*ILENG
C     J40:---
      J50=J40
C     J50:---
      J90=J50
C     J90:KCP(MAXGRD)
      J100=J90 +MAXGRD
C     J100:COEF(3*ILENG)
      J110=J100+3*ILENG
C     J110:DTRAI_(4*ILENG)
      J120=J110+4*ILENG
C     J120:DRHOI_(4*MAXGRD)
      J130=J120+4*MAXGRD
C     J130:EX0(ILENG)
      J140=J130+ILENG
C     J140:EX(18*ILENG)
      J150=J140+18*ILENG
C     J150:EC0(ILENG)
      J160=J150+ILENG
C     J160:EC(55*ILENG)
      J170=J160+55*ILENG
      J180=J170
C     J180:AOMAX(L1)
      J190=J180+L1
C     J190:IAO(L1)
      J200=J190+L1
C     J200:IPRGRD(2*L3)
      LAST=J200+2*L3
      ELSE
      LAST=LASTI
      ENDIF
      IDUM=LAST
      NEED=LAST-I10
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 800
C
      IF (NDFTFG.EQ.1) THEN
C
C     ----- CONSTRUCT DENSITY
C
      CALL TDDFTSET(RHOSET,X(J10),X(J20),X(J30),V,
     *              X(J120),X(J180),X(J190),ILENG,NPTGRD,L1)
C        ROUTINE TDDFTSET(RHOSET,XYZGRD,XYZWGT,XYZDCH,V,
C    *                    DRHOI_,AOMAX,GMO,ILENG,NPTGRD,L1)
      ENDIF
C
C     --- GET M(IA,JB) AS (EA-EI)
C
      CALL XECALC(X(I20),E,NOCC,LX,L7)
C
C     --- GET M^(-1)
C
      CALL PCGMIN(X(I20),X(I30),L7)
C     SUBROUTINE PCGMIN(XM,XM_IN,L7)
C
C     --- GET INITIAL X_T(L7) (= ZERO)
C
      CALL VCLR(X(I40),1,L7)
C     TRANSFORM INTO AO -BASIS
      CALL IATOGENTD(X(I40),X(K10),LX,L7,NOCC,1,1)
C     SUBROUTINE IATOGENTD(PV,AV,LX,L7,NOCC,IVEC,MXVEC)
      CALL DGEMM('N','N',L1,LX,LX,ONE,V,L1,X(K10),LX,ZERO
     * ,X(I100),L1)
      CALL DGEMM('N','T',L1,L1,LX,ONE,X(I100),L1,V,L1,ZERO
     * ,X(K10),L1)
C
C     --- GET INITIAL A_PK(L7)
C
C ********************************************************************
C
C     ---- TWO ELECTRON PART
C
      CALL VCLR(X(K20),1,L2)
      CALL TD2E(X(K10),X(K20),X(IDUM),X(K40),X(K50),L1,ONE,0,1,1)
C
C     -- XC
C
      IF(NDFTFG.EQ.1) THEN
         CALL VCLR(X(K30),1,L2)
         CALL TDFXCP(X(K30),RHO,X(J10),X(J20),X(J30),X(K10),X(IDUM),
     *               X(J120),X(J110),X(J100),X(J140),X(J160),X(J130),
     *               X(J150),X(J180),X(J200),X(J190),ILENG,NPTGRD,
     *               L1,L2,LX,2,2)
         IF(GOPARR.AND.IBTYP.EQ.1) THEN
            CALL DDI_GSUMF(2338,X(K30), L2)
            CALL DDI_GSUMF(2339,RHO,1)
         ENDIF
         CALL DAXPY(L2,2.0D+00,X(K30),1,X(K20),1)
      ENDIF
C
C     -- PRODUCT A_P
C
      CALL CANTOGEN(X(K20),X(K10),L1)
      CALL MNTOIA(X(K10),X(I50),V,X(I100),L1,LX,L7,NOCC,1)
      CALL ESUM(E,X(I50),X(I40),LX,L7,NOCC,1)
C
C     --- GET R,BETA,PK INITIAL
C     USE R(Z)
      CALL PCGRBPINI(R,X(I30),X(I50),X(I70),BETA,X(I80),L7,ER)
C     SUBROUTINE PCGRBPINI(D,XM_IN,A_PK,R,BETA,PK,L7,ER)
      IF(MASWRK) WRITE(IW,1111) ER,THRE
      IF (ER.LE.THRE) THEN
       WRITE(IW,*)'ACCIDENTAL CONVERGENCE!'
       ITER=0
       GOTO 100
       ELSE
       CONTINUE
      ENDIF
C
C     --- START PCG LOOP ---
C
      ITER=0
   10 CONTINUE
      ITER=ITER+1
      IF(ITER.GT.MXITER) THEN
       IF(MASWRK) THEN
       WRITE(6,*)'    XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX   '
       WRITE(6,*)'      Z-VECTOR FAILED TO CONVERGE'
       WRITE(6,*)'    XXXXXXXXXXXX STOP XXXXXXXXXXXXX   '
       WRITE(6,*)' ITER/MXITER=',ITER,'/',MXITER
       ENDIF
       CALL ABRT
       STOP
      ENDIF
C
C     -- TRANSFORM TRIAL VECTOR INTO AO BASIS
C
      CALL IATOGENTD(X(I80),X(K10),LX,L7,NOCC,1,1)
C     SUBROUTINE IATOGENTD(PV,AV,LX,L7,NOCC,IVEC,MXVEC)
      CALL DGEMM('N','N',L1,LX,LX,ONE,V,L1,X(K10),LX,ZERO
     * ,X(I100),L1)
      CALL DGEMM('N','T',L1,L1,LX,ONE,X(I100),L1,V,L1,ZERO
     * ,X(K10),L1)
C
C     -- GET A_PK
C
C ********************************************************************
C
C     ---- TWO ELECTRON PART
C
      CALL VCLR(X(K20),1,L2)
      CALL TD2E(X(K10),X(K20),X(IDUM),X(K40),X(K50),L1,ONE,0,1,1)
C
C     -- XC
C
      IF(NDFTFG.EQ.1) THEN
         CALL VCLR(X(K30),1,L2)
         CALL TDFXCP(X(K30),RHO,X(J10),X(J20),X(J30),X(K10),X(IDUM),
     *               X(J120),X(J110),X(J100),X(J140),X(J160),X(J130),
     *               X(J150),X(J180),X(J200),X(J190),ILENG,NPTGRD,
     *               L1,L2,LX,2,2)
         IF(GOPARR.AND.IBTYP.EQ.1) THEN
            CALL DDI_GSUMF(2340,X(K30),L2)
            CALL DDI_GSUMF(2341,RHO,1)
         ENDIF
         CALL DAXPY(L2,2.0D+00,X(K30),1,X(K20),1)
      ENDIF
C
C     -- PRODUCT A_P
C
       CALL CANTOGEN(X(K20),X(K10),L1)
       CALL MNTOIA(X(K10),X(I50),V,X(I100),L1,LX,L7,NOCC,1)
       CALL ESUM(E,X(I50),X(I80),LX,L7,NOCC,1)
C ********************************************************************
C
C     -- GET ALPHA
C
      SUM=DDOT(L7,X(I80),1,X(I50),1)
      ALPHA=1.0D+00/SUM
C
C     -- UPDATE VECTOR X : X=ALPHA*PK
C
      CALL DAXPY(L7,ALPHA,X(I80),1,X(I40),1)
C
C     -- GET ERROR VECTOR R(L7) & CHECK CONVERGENCE
C
      CALL DAXPY(L7,-ALPHA,X(I50),1,X(I70),1)
      ER=DDOT(L7,X(I70),1,X(I70),1)
      IF(MASWRK) WRITE(IW,1112) ITER,ER,THRE
      IF (ER.LE.THRE) GOTO 100
C
C     -- GET BETA
C
      CALL PCGB(BETA,X(I70),X(I30),L7)
C     SUBROUTINE PCGB(BETA,R,XM_IN,L7)
C
C     -- UPDATE PK(L7)
C
      CALL PCGPUP(X(I80),BETA,X(I30),X(I70),L7)
C     SUBROUTINE PCGPUP(PK,BETA,XM_IN,R,L7)
      GOTO 10
C
  100 CONTINUE
      IF(MASWRK) THEN
      WRITE(IW,*)'   -------------------------'
      WRITE(IW,*)'       Z-VECTOR CONVERGED          '
      WRITE(IW,*)'   -------------------------'
      ENDIF
C
C     --- CONVERGED
C     NOW, X_T(L7) EQUAL TO ZIA(L7)
      CALL DCOPY(L7,X(I40),1,R,1)
C
C     ----- CONVERT Z_IA TO Z_MN
C
      CALL IATOGENTD(X(I40),X(I100),LX,L7,NOCC,1,1)
C     SUBROUTINE IATOGENTD(PV,AV,LX,L7,NOCC,IVEC,MXVEC)
      CALL DGEMM('N','N',L1,LX,LX,ONE,V,L1,X(I100),LX,ZERO
     * ,X(K10),L1)
      CALL DGEMM('N','T',L1,L1,LX,ONE,X(K10),L1,V,L1,ZERO
     * ,X(I100),L1)
C     -- PA=TA+ZA
      CALL DAXPY(L1*L1,1.0D+00,X(I100),1,TA,1)
      CALL DAWRIT(IDAF,IODA,TA,L3,IRECTD+4,0)
C---  IF(MASWRK) WRITE(IW,*) 'PA WRIT DONE'
C
  800 CONTINUE
      CALL RETFM(NEED)
      RETURN
 1111 FORMAT(' INITIAL ERROR =',3X,1P,E10.3,1X,'/',1P,E10.3)
 1112 FORMAT(' ITER#',I2,' ERROR =',3X,1P,E10.3,1X,'/',1P,E10.3)
      END
C*MODULE TDDGRD  *DECK XECALC
      SUBROUTINE XECALC(XE,E,NOCC,LX,L7)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION XE(L7)
      DIMENSION E(LX)
      DO 10 I=1,NOCC
      DO 10 J=NOCC+1,LX
      JJ=J-NOCC
      IJ=(JJ-1)*NOCC+I
      XE(IJ)=E(J)-E(I)
   10 CONTINUE
      RETURN
      END
C*MODULE TDDGRD  *DECK PCGB
      SUBROUTINE PCGB(BETA,R,XM_IN,NDIM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION R(*)
      DIMENSION XM_IN(*)
      DUMI=0.0D+00
      DO 10 I=1,NDIM
      DUM=R(I)*XM_IN(I)*R(I)
   10 DUMI=DUMI+DUM
      BETA=1.0D+00/DUMI
      RETURN
      END
C*MODULE TDDGRD  *DECK PCGPUP
      SUBROUTINE PCGPUP(PK,BETA,XM_IN,R,NDIM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION PK(*)
      DIMENSION XM_IN(*)
      DIMENSION R(*)
      DO 10 I=1,NDIM
   10 PK(I)=PK(I)+BETA*XM_IN(I)*R(I)
      RETURN
      END
C*MODULE TDDGRD  *DECK PCGMIN
      SUBROUTINE PCGMIN(XM,XM_IN,NDIM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION XM(*)
      DIMENSION XM_IN(*)
      DO 10 I=1,NDIM
   10 XM_IN(I)=1.0D+00/XM(I)
      RETURN
      END
C*MODULE TDDGRD  *DECK PCGRBPINI
      SUBROUTINE PCGRBPINI(D,XM_IN,A_PK,R,BETA,PK,NDIM,ER)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(*)
      DIMENSION XM_IN(*)
      DIMENSION A_PK(*)
      DIMENSION R(*)
      DIMENSION PK(*)
C
C     -- R INI AND R NORM
C
      RNORM=0.0D+00
      DO 10 I=1,NDIM
      R(I)=D(I)-A_PK(I)
   10 RNORM=RNORM+R(I)**2
      ER=RNORM
C
C     -- BETA INI
C
      DUMI=0.0D+00
      DO 20 I=1,NDIM
      DUM=R(I)*XM_IN(I)*R(I)
   20 DUMI=DUMI+DUM
      DUM=1.0D+00/DUMI
      BETA=DUM
C
C     -- PK INI
C
      DO 30 I=1,NDIM
   30 PK(I)=BETA*XM_IN(I)*R(I)
C
      RETURN
      END
C*MODULE TDDGRD  *DECK WMNG
      SUBROUTINE WMNG(QAI,GXC,P,Z,V1,V2,V,E,NOCC,L0,L1,L2,L3,LX,
     *                LABSIZ,NEED)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,TRIPLET,SG1,SG1T
C
      DIMENSION QAI(*),GXC(*),P(*),Z(*),V1(*),V2(*),V(*),E(*)
C
      PARAMETER (MXRT=100)
      PARAMETER (MXGRID=10)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /FMCOM / X(1)
      COMMON /INFGRD/ RHOMIN,ILENG,MAXGRD
      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 /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
      CALL VALFM(LOADFM)
      CALL GOTFM(NGOTMX)
      I10 = LOADFM + 1
C     I10:HPP_AO(L2)
      I15 = I10+L2
C     I15:PA(L3) OR SCRATCH
      I20=I15+L3
C     I20:XX(L0)
      I30=I20+L0
C     I30:IX(L0*LABSIZ)
      I40=I30+L0*LABSIZ
C     I40:--
      I50=I40
C     I50:--
      I60=I50
C     I60:HPPIJ(NOCC*NOCC)
      I70=I60+NOCC*NOCC
C     I70:FXCP(L2)
      I80=I70+L2
C     I80:GIJ(NOCC,NOCC)
      I90=I80+NOCC*NOCC
C     I90:W(LX,LX)
      I100=I90+LX*LX
C     I100:--
      I110=I100
C     I110:W_AO(L3) OR SCRATCH
      LASTI=I110+L3
C
C     ----- MEMORY FOR DFT GRID
C
      IF(NDFTFG.EQ.1) THEN
      J10=LASTI+1
C     J10:XYZGRD(MAXGRD*3)
      J20=J10+MAXGRD*3
C     J20:XYZWGT(MAXGRD)
      J30=J20+MAXGRD
C     J30:XYZDCH
      J40=J30+4*L1*ILENG
C     J40:---
      J50=J40
C     J50:---
      J90=J50
C     J90:KCP(MAXGRD)
      J100=J90 +MAXGRD
C     J100:COEF(3*ILENG)
      J110=J100+3*ILENG
C     J110:DTRAI_(4*ILENG)
      J120=J110+4*ILENG
C     J120:DRHOI_(4*MAXGRD)
      J130=J120+4*MAXGRD
C     J130:EX0(ILENG)
      J140=J130+ILENG
C     J140:EX(18*ILENG)
      J150=J140+18*ILENG
C     J150:EC0(ILENG)
      J160=J150+ILENG
C     J160:EC(55*ILENG)
      J170=J160+55*ILENG
      J180=J170
C     J180:AOMAX(L1)
      J190=J180+L1
C     J190:IAO(L1)
      J200=J190+L1
C     J200:IPRGRD(2*L3)
      LAST=J200+2*L3
      ELSE
      LAST=LASTI
      ENDIF
      IDUM=LAST
      NEED=LAST-I10
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C     ----- FPP(L2) @ X(I10)
C
C     ---- TWO ELECTRON PART
C
      CALL VCLR(X(I10),1,L2)
      CALL TD2E(P,X(I10),X(IDUM),X(I20),X(I30),L1,ONE,0,1,1)
C
C     -- XC
C
      IF(NDFTFG.EQ.1) THEN
         CALL TDDFTSET(RHOSET,X(J10),X(J20),X(J30),V,X(J120),
     *                 X(J180),X(J190),ILENG,NPTGRD,L1)
         CALL TDFXCP(X(I70),RHO,X(J10),X(J20),X(J30),P,X(IDUM),
     *               X(J120),X(J110),X(J100),X(J140),X(J160),
     *               X(J130),X(J150),X(J180),X(J200),X(J190),
     *               ILENG,NPTGRD,L1,L2,LX,2,2)
         IF(GOPARR.AND.IBTYP.EQ.1) THEN
            CALL DDI_GSUMF(2342,X(I70), L2)
            CALL DDI_GSUMF(2343,RHO,1)
         ENDIF
         CALL DAXPY(L2,2.0D+00,X(I70),1,X(I10),1)
      ENDIF
C
C     ---- HPPAO TO HPPIJ @ X(I60)
C
      CALL CANTOGEN(X(I10),X(I15),L1)
      CALL DGEMM('N','N',L1,NOCC,L1,ONE,X(I15),L1,V,L1,ZERO
     * ,X(I110),L1)
      CALL DGEMM('T','N',NOCC,NOCC,L1,ONE,V,L1,X(I110),L1,ZERO
     * ,X(I60),NOCC)
C     ---- GAO TO GIJ
      IF(NDFTFG.EQ.1) THEN
       CALL CANTOGEN(GXC,X(I15),L1)
       CALL DGEMM('N','N',L1,NOCC,L1,ONE,X(I15),L1,V,L1,ZERO
     * ,X(I110),L1)
       CALL DGEMM('T','N',NOCC,NOCC,L1,ONE,V,L1,X(I110),L1,ZERO
     * ,X(I80),NOCC)
       ELSE
       CALL VCLR(X(I80),1,NOCC*NOCC)
       ENDIF
C
C       ---- CALCULATE W_MO
C
      CALL WCALC(X(I90),ESTATE(NTHST),E,V1,V2,Z,LX,NOCC,NOCC,
     * QAI,X(I80),X(I60))
C
C     ---- W_MO TO W_AO AND SAVE W_AO
C
      CALL DGEMM('N','N',L1,LX,LX,ONE,V,L1,X(I90),LX,ZERO
     * ,X(I15),L1)
      CALL DGEMM('N','T',L1,L1,LX,ONE,X(I15),L1,V,L1,ZERO
     * ,X(I110),L1)
      CALL DAWRIT(IDAF,IODA,X(I110),L3,IRECTD+5,0)
C     IF(MASWRK) WRITE(IW,*)'WA WRIT DONE'
C
  800 CONTINUE
      CALL RETFM(NEED)
C
      RETURN
      END
C*MODULE TDDGRD  *DECK WCALC
      SUBROUTINE WCALC(W,EE,E,V1,V2,Z,LX,NOCCG,NOCC,
     *                 XHX,GIJ,HPPIJ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION W(LX,*),E(*),V1(*),V2(*),Z(*),XHX(*),GIJ(NOCCG,*),
     *          HPPIJ(NOCCG,*)
C
      LOGICAL SG1
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
C
      DATA ZERO/0.0D+00/
C
C     ---- W_IJ
C
      CALL VCLR(W,1,LX*LX)
      DO I=1,NOCC
      DO J=1,I
C     /--FOR DUM1
      DUM1=ZERO
      DO K=NOCC+1,LX
      KK=K-NOCC
      IIA=(KK-1)*NOCC+I
      JA=(KK-1)*NOCC+J
        DUM=V1(IIA)*V2(JA)+V2(IIA)*V1(JA)
      DUM1=DUM1+DUM
      END DO
      DUM1=EE*DUM1
C     /--FOR DUM2 ----
      DUM2=ZERO
      DO K=NOCC+1,LX
      KK=K-NOCC
      IIA=(KK-1)*NOCC+I
      JA=(KK-1)*NOCC+J
      DUM=V1(IIA)*V1(JA)+V2(IIA)*V2(JA)
      DUM2=DUM2+E(K)*DUM
      END DO
C     /-- FOR DUM3 -----
        DUM3= HPPIJ(I,J)
        W(I,J)=DUM1-DUM2+DUM3
C     --- ADD GIJ TO WIJ (DFT)
        IF(NDFTFG.EQ.1) W(I,J)=2.0D+00*GIJ(I,J)+W(I,J)
      END DO
      END DO
C
C     ---- W_AB
C
      DO I=NOCC+1,LX
      DO J=NOCC+1,I
      II=I-NOCC
      JJ=J-NOCC
C     ---FOR DUM1
        DUM1=ZERO
        DO K=1,NOCC
      IIA=(II-1)*NOCC+K
      IB=(JJ-1)*NOCC+K
        DUM=V1(IIA)*V2(IB)+V2(IIA)*V1(IB)
      DUM1=DUM1+DUM
      END DO
      DUM2=ZERO
      DO K=1,NOCC
      IIA=(II-1)*NOCC+K
      IB=(JJ-1)*NOCC+K
      DUM=E(K)*(V1(IIA)*V1(IB)+V2(IIA)*V2(IB))
      DUM2=DUM2+DUM
      END DO
      W(I,J)=EE*DUM1+DUM2
      END DO
      END DO
C
C     ---- W_IA
C
      DO I=1,NOCC
      DO J=NOCC+1,LX
      JJ=J-NOCC
      IJX=(JJ-1)*NOCC+I
      W(J,I)=XHX(IJX)+E(I)*Z(IJX)
      END DO
      END DO
C
C     ---- SCALE HALF TO DIAGONAL ELEMENTS
C
      DO I=1,LX
      W(I,I)=0.5D+00*W(I,I)
      END DO
C
      CALL DSCAL(LX*LX,-1.0D+00,W,1)
C
      RETURN
      END
C*MODULE TDDGRD  *DECK DENTD
      SUBROUTINE DENTD(DA,P,L1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL TRIPLET,SG1T
      PARAMETER (MXRT=100)
      DIMENSION DA(*),P(L1,L1)
C
      COMMON /IJPAIR/ IA(8192)
      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)
C
C     -- READ P
      CALL DAREAD(IDAF,IODA,P,L1*L1,IRECTD+4,0)
C     -- TOTAL DENSITY
      DO 10 M=1,L1
      DO 10 N=1,M
      MN=IA(M)+N
      PDUM=P(M,N)+P(N,M)
   10 DA(MN)=DA(MN)+PDUM
      RETURN
      END
C*MODULE TDDGRD  *DECK TDDENS
      SUBROUTINE TDDENS(EPS,W,L1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL TRIPLET,SG1T
      DIMENSION EPS(*),W(L1,*)
      PARAMETER (MXRT=100)
C
      COMMON /IJPAIR/ IA(8192)
      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)
C
C     -- READ W
      CALL DAREAD(IDAF,IODA,W,L1*L1,IRECTD+5,0)
      DO 10 M=1,L1
      DO 10 N=1,M
      MN=IA(M)+N
   10 EPS(MN)=W(M,N)+W(N,M)+EPS(MN)
      RETURN
      END
C*MODULE TDDGRD  *DECK TDGAMMA
      SUBROUTINE TDGAMMA(M,N,K,L,D,P,V1,V2,SD,L1,CEXC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL LCFLAG,LRINT
      COMMON /IJPAIR/ IA(8192)
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      PARAMETER(PT5=0.5D+00,FOUR=4.0D+00)
      DIMENSION V1(L1,*),V2(L1,*),P(L1,*),D(*)
C     -- SET INDICES
      KL=IA(MAX0(K,L))+MIN0(K,L)
C     MN=IA(MAX0(M,N))+MIN0(M,N)
C     ML=IA(MAX0(M,L))+MIN0(M,L)
C     MK=IA(MAX0(M,K))+MIN0(M,K)
      KN=IA(MAX0(K,N))+MIN0(K,N)
      LN=IA(MAX0(L,N))+MIN0(L,N)
C     -- P=T+Z CONTRIBUTION
      DUM1C=FOUR*P(M,N)*PT5*D(KL)
      DUM1E=-CEXC*(P(M,L)*PT5*D(KN)+P(K,M)*PT5*D(LN))
C     -- OTHER MO-INTEGRAL CONTRIBUTION
      DUM2C=FOUR*V1(M,N)*V1(K,L)
      DUM2=V1(K,N)*V1(M,L)
     * +V1(K,M)*V1(L,N)-V2(K,N)*V2(M,L)+V2(K,M)*V2(L,N)
      DUM2E=-CEXC*DUM2
C     -- TOTAL
      SD= DUM1C+DUM1E + DUM2C+DUM2E
      IF(LCFLAG.AND.LRINT) SD=DUM1E+DUM2E
      IF(LCFLAG.AND.(.NOT.LRINT)) SD=DUM1C+DUM2C
      SD= FOUR*SD
C     SD=TWO*D(MN)*D(KL)-D(MK)*D(LN)
C
      RETURN
      END
C*MODULE TDDGRD  *DECK JKDRPA
      SUBROUTINE JKDRPA(DA,DB,P,V1,V2,L3,L2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION DA(*),DB(*),P(*),V1(*),V2(*)
C
      LOGICAL TRIPLET,SG1T
C
      PARAMETER (MXRT=100)
C
      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,IJK,IPK,IDAF,NAV,IODA(950)
C
      CALL DAREAD(IDAF,IODA,DA,L2,308,0)
      CALL DAREAD(IDAF,IODA,DB,L2,308,0)
      CALL DAREAD(IDAF,IODA, P,L3,IRECTD+4,0)
      CALL DAREAD(IDAF,IODA,V1,L3,IRECTD+2,0)
      CALL DAREAD(IDAF,IODA,V2,L3,IRECTD+3,0)
      RETURN
      END
C*MODULE TDDGRD  *DECK TDDER
      SUBROUTINE TDDER
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXRT=100)
C
      LOGICAL TRIPLET,SG1T,DEBUG,GOPARR,DSKWRK,MASWRK
C
      COMMON /FMCOM / X(1)
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /INFGRD/RHOMIN,ILENG,MAXGRD
      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 /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK,DBUGME/8HCHECK   ,8HTDDER   /
C
C         COMPUTE GRID POINT CONTRIBUTIONS TO TD-DFT NUCLEAR GRADIENT
C
      DEBUG = EXETYP.EQ.DBUGME
      IF(MASWRK  .AND.  DEBUG) THEN
         WRITE(IW,*) 'GRADIENT VECTOR ON ENTRY TO TDDER'
         CALL EGOUT(DE,NAT)
      END IF
C
C     --- DEFINE MAXGRD
C
      CALL MEMGRD
C
C     --- MEMORY
C
      L1= NUM
      L2=(NUM*(NUM+1))/2
      L3= NUM* NUM
      LX= NQMT
      CALL VALFM(LOADFM)
      I10=LOADFM+1
C     I10:V(L1*LX) OR SCRATCH
      I20=I10+L1*L1
C     I20:PA(L3)
      I30=I20+L3
C     I30:V1A(L3)
      J10=I30+L3
C     J10:XYZGRD(MAXGRD*3)
      J20=J10+MAXGRD*3
C     J20:XYZWGT(MAXGRD)
      J30=J20+MAXGRD
C     J30:XYZDDC(10,NUM,ILENG)
      J40=J30+10*L1*ILENG
C     J40:D(L2)
      J50=J40+L2
C     J50:KCP(MAXGRD)
      J60=J50+MAXGRD
C     J60:GRDFUN(MAXGRD)
      J70=J60+MAXGRD
C     J70:DRHOI_(4*MAXGRD)
      J80=J70+4*MAXGRD
C     J80:EX0(ILENG)
      J90=J80+ILENG
C     J90:EX(18*ILENG)
      J100=J90+18*ILENG
C     J100:EC0(ILENG)
      J110=J100+ILENG
C     J110:EC(55*ILENG)
      J120=J110+55*ILENG
C     J120:AOMAX(L1)
      J130=J120+L1
C     J130:IPRGRD(2,L3)
      J140=J130+2*L3
C     J140:IAO(L1)
      J150=J140+L1
C     J150:DDUM(ILENG,4)
      J160=J150+4*ILENG
C     J160:PDUM(ILENG,4)
      J170=J160+4*ILENG
C     J170:VDUM(ILENG,4)
      J180=J170+4*ILENG
C     J180:DTRAI_(4*ILENG)
      LAST=J180+4*ILENG
      NEED=LAST-I10
      CALL GETFM(NEED)
C
C        ADDITIONAL MEMORY IS ALLOCATED INSIDE TDWGT, EXECUTE IT TOO.
C
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C     -- GRID READ AND DENSITY SET
C
      CALL DAREAD(IDAF,IODA,X(I10),L1*LX,15,0)
      CALL TDDFTSET(RHOSET,X(J10),X(J20),X(J30),X(I10),
     *              X(J70),X(J120),X(J130),ILENG,NPTGRD,L1)
C
C     ----- CALCULATE XC GRADIENT     BY USING PA, V1A, DA
C
      CALL DAREAD(IDAF,IODA,X(I20),L1*L1,IRECTD+4,0)
      CALL DAREAD(IDAF,IODA,X(I30),L1*L1,IRECTD+2,0)
      CALL DAREAD(IDAF,IODA,X(J40),L2,308,0)
      CALL DSCAL(L2,0.5D+00,X(J40),1)
      CALL CANTOGEN(X(J40),X(I10),L1)
      CALL VFEXC(X(I10),X(I20),X(I30),X(J10),X(J20),X(J30),X(J70),
     *           X(J180),X(J80),X(J90),X(J100),X(J110),X(J60),X(J50),
     *           X(J120),X(J150),X(J160),X(J170),X(J140),X(J130),
     *           NPTGRD,L1,ILENG)
      IF(MASWRK  .AND.  DEBUG) THEN
         WRITE(IW,*) 'GRADIENT VECTOR AFTER VFEXC IN TDDER'
         CALL EGOUT(DE,NAT)
      END IF
C
C     ----- WEIGHT DERIVATIVE CORRECTION
C
  800 CONTINUE
      NEEDDER = NEED
      CALL TDWGT(X(J60),X(J70),NPTGRD,NEEDDER)
      IF(MASWRK  .AND.  DEBUG) THEN
         WRITE(IW,*) 'GRADIENT VECTOR AFTER TDWGT IN TDDER'
         CALL EGOUT(DE,NAT)
      END IF
C
      CALL RETFM(NEED)
      CALL TIMIT(1)
      RETURN
      END
C*MODULE TDDGRD  *DECK VFEXC
      SUBROUTINE VFEXC(D,PA,V1A,XYZGRD,XYZWGT,XYZDDC,DRHOI_,
     *                 DTRAI_,EX0,EX,EC0,EC,GRDFUN,KCP,AOMAX,
     *                 DDUM,PDUM,VDUM,IAO,IPRGRD,NPTGRD,L1,ILENG)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION D(*),PA(*),V1A(*),
     *          XYZGRD(3,NPTGRD),XYZWGT(NPTGRD),XYZDDC(*),
     *          DRHOI_(*),DTRAI_(*),EX0(*),EX(*),EC0(*),EC(*),
     *          GRDFUN(NPTGRD),KCP(NPTGRD),AOMAX(*),
     *          DDUM(*),PDUM(*),VDUM(*),IAO(*),IPRGRD(*)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT,DBUG,SG1
C
      PARAMETER (MXATM=2000)
      PARAMETER (MXGRID=10,MXGRIDTYP=10)
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 /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
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
      DBUG =.FALSE.
      NPT=NRAD*MAXANG
      RCUTOFF=CONVHF/(NPT*NAT)
C
      CALL GRDRD(22,XYZGRD,XYZWGT,KCP,NPTGRD,1)
      CALL VCLR(GRDFUN,1,NPTGRD)
      RHO=0.0D+00
      IST=0
      IEND=0
      IAREA=0
C
C     LOAD BALANCING TERMS
      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 DO LOOP FOR PARALLELIZATION
C
      DO 10 ILG=1,NLENG
C   10 CONTINUE
C
C     DYNAMIC AND STATIC LOAD BALANCING
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
C      IST=IEND+1
      IST=(ILG-1)*ILENG+1
      IEND=IST+ILENG-1
      IAREA=IAREA+1
      IF(IEND.GT.NPTGRD) IEND=NPTGRD
C
C     ---- FUNCTIONAL ENERGY AND DERIVATIVES
C          EX(ILENG,18) & EC(ILENG,55)
      CALL DFXCSET(EX0,EX,EC0,EC,DRHOI_,XYZWGT,1.0D-12,EXC,ECOR,
     * ILENG,NPTGRD,1,3,IST,IEND)
C
C     ---- GRADIENT DENSIY:DDUM,PDUM,VDUM
C
      CALL GRDGAO(IST,IEND,XYZGRD,XYZDDC,AOMAX,L1)
      THRE=1.0D-06
      CALL AOSERCH(THRE,AOMAX,IAO,L1,IAON)
C     PA
      THRE=1.0D-11
      CALL SELECTPR(THRE,PA,AOMAX,IPRGRD,IAO,IMAX,L1,IAON)
      CALL TRADEN2(RCUTOFF,DRHOI_,DTRAI_,PA,XYZDDC,IPRGRD,
     * IST,IEND,NPTGRD,IMAX,L1,ILENG)
      CALL GSVCOEF(DDUM,PDUM,VDUM,DRHOI_,DTRAI_,EX,EC,EX0,EC0
     * ,GRDFUN,RCUTOFF,NPTGRD,ILENG,IST,IEND,1)
C     V1A
      CALL SELECTPR(THRE,V1A,AOMAX,IPRGRD,IAO,IMAX,L1,IAON)
      CALL TRADEN2(RCUTOFF,DRHOI_,DTRAI_,V1A,XYZDDC,IPRGRD,
     * IST,IEND,NPTGRD,IMAX,L1,ILENG)
      CALL FCOEF(DDUM,VDUM,DRHOI_,DTRAI_,EX,EC
     * ,GRDFUN,RCUTOFF,NPTGRD,ILENG,IST,IEND,1)
C
C     ---- V & F GRADIENT
C
      CALL GRDGGAO(IST,IEND,XYZGRD,XYZDDC,AOMAX,L1)
      CALL XCGRD(RHO,DRHOI_,DDUM,PDUM,VDUM,D,D,PA,V1A,
     * XYZWGT,XYZDDC,RCUTOFF,NPTGRD,KCP,ILENG,IST,IEND,L1,1)
C
C      IF(IEND.EQ.NPTGRD) THEN
C      GOTO 20
C      ELSE
C      GOTO 10
C      ENDIF
C   20 CONTINUE
C
 10   CONTINUE
C
C     RESET FOR DYNAMIC LOAD BALANCING
C
      IF(GOPARR.AND.NXT) CALL DDI_DLBRESET
C
      IF(DBUG) WRITE(6,*)'ELECTRONS IN XCGRAD=',RHO
      RETURN
      END
C*MODULE TDDGRD  *DECK TDWGT
      SUBROUTINE TDWGT(GRDFUN,DRHOI_,NPTGRD,NEEDDER)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      LOGICAL GOPARR,MASWRK,DSKWRK
C
      DIMENSION GRDFUN(*),DRHOI_(*)
C
      COMMON /FMCOM / X(1)
      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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      CALL VALFM(LOADFM)
      CALL GRDDFT(L2,NGRDMEM)
      IGRDDFT = LOADFM  + 1
      IUVEC   = IGRDDFT + NGRDMEM
      IDWTINT = IUVEC   + NAT*3
      IDWTTOT = IDWTINT + 3*NAT*NAT
      IDTOTWT = IDWTTOT + 3*NAT
      LAST    = IDTOTWT + 3*NAT
      NEED = LAST-LOADFM-1
      IF(MASWRK) WRITE(IW,9010) NEED+NEEDDER
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C     -- DFTSET
C
      CALL DFTSET(X(IGRDDFT),1,.FALSE.)
C
C     --- WEIGHT DERIVATIVE CORRECTION
C     PARALLELIZATION IS TURNED OFF WITHIN WGTGRAD, DUE TO DERWGT
C
      CALL WGTGRAD(X(IGRDDFT),X(IUVEC),X(IDWTINT),X(IDWTTOT),X(IDTOTWT),
     *             TOTELE,DRHOI_,GRDFUN,NPTGRD,IGPT)
C
      IF(IGPT.NE.NPTGRD) THEN
         IF(MASWRK) THEN
            WRITE(IW,*) 'SOMETHING IS WRONG.. IGPT.NE.NPTGRD'
            WRITE(IW,*) 'IGPT/NPTGRD=',IGPT,'/',NPTGRD
            WRITE(IW,*) 'NUMBER OF ELECTRONS=',TOTELE
         END IF
         CALL ABRT
      ENDIF
C
  800 CONTINUE
      CALL RETFM(NEED)
      RETURN
C
 9010 FORMAT(1X,'MEMORY FOR GRID POINT CONTRIBUTIONS TO THE',
     *          ' TD-DFT GRADIENT=',I10,' WORDS.')
      END
C*MODULE TDDGRD  *DECK XCGRD
      SUBROUTINE XCGRD(RHO,DRHOI_,DDUM,PDUM,VDUM,DA,DB,PA,V1A,
     *                 XYZWGT,XYZDDC,RHOMIN,NPTGRD,KCP,
     *                 ILENG,IST,IEND,L1,IRU)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXGTOT=20000,MXATM=2000)
C
      DIMENSION DA(L1,*),DB(L1,*)
      DIMENSION PA(L1,L1,IRU)
      DIMENSION V1A(L1,L1,IRU)
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION DDUM(ILENG,4,IRU)
      DIMENSION PDUM(ILENG,4,IRU)
      DIMENSION VDUM(ILENG,4,IRU)
      DIMENSION XYZWGT(*)
      DIMENSION XYZDDC(10,L1,*)
      DIMENSION KCP(*)
      DIMENSION DUMA(4,2),DEMA(4,2),DOMA(4,2)
C
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
C
C     -- GRID LOOP
      DO IPT=IST,IEND
      IIPT=IPT-IST+1
      RHOA = DRHOI_(IPT,1,1)
      RHOB = DRHOI_(IPT,1,IRU)
      RHOT=RHOA+RHOB
      IF(RHOT.GT.RHOMIN) THEN
C
C     --- GRADIENT OF THE THREE DIFFERENT DENSITIES
C     -- D DENSITY GRADIENT
      DUMA(1,1) = DDUM(IIPT,1,1)
      DUMA(2,1) = DDUM(IIPT,2,1)
      DUMA(3,1) = DDUM(IIPT,3,1)
      DUMA(4,1) = DDUM(IIPT,4,1)
      DUMA(1,2) = DDUM(IIPT,1,IRU)
      DUMA(2,2) = DDUM(IIPT,2,IRU)
      DUMA(3,2) = DDUM(IIPT,3,IRU)
      DUMA(4,2) = DDUM(IIPT,4,IRU)
C     -- P DENSITY GRADIENT
      DEMA(1,1) = PDUM(IIPT,1,1)
      DEMA(2,1) = PDUM(IIPT,2,1)
      DEMA(3,1) = PDUM(IIPT,3,1)
      DEMA(4,1) = PDUM(IIPT,4,1)
      DEMA(1,2) = PDUM(IIPT,1,IRU)
      DEMA(2,2) = PDUM(IIPT,2,IRU)
      DEMA(3,2) = PDUM(IIPT,3,IRU)
      DEMA(4,2) = PDUM(IIPT,4,IRU)
C     -- V1 DENSITY GRADIENT
      DOMA(1,1) = VDUM(IIPT,1,1)
      DOMA(2,1) = VDUM(IIPT,2,1)
      DOMA(3,1) = VDUM(IIPT,3,1)
      DOMA(4,1) = VDUM(IIPT,4,1)
      DOMA(1,2) = VDUM(IIPT,1,IRU)
      DOMA(2,2) = VDUM(IIPT,2,IRU)
      DOMA(3,2) = VDUM(IIPT,3,IRU)
      DOMA(4,2) = VDUM(IIPT,4,IRU)
C     CALL VCLR(DUMA,1,4*2)
C     CALL VCLR(DEMA,1,4*2)
C     CALL VCLR(DOMA,1,4*2)
C
C     --- CALCULATE DENSITY GRADIENT
C
      KAT=KCP(IPT)
      DO 100 ISHELL=1,NSHELL
       IAT =KATOM(ISHELL)
       IF(IAT.EQ.KAT) GOTO 100
       MINI=KMIN (ISHELL)
       MAXI=KMAX (ISHELL)
       LOCI=KLOC (ISHELL)-MINI
       DO IBFN=MINI,MAXI
       ICHI=LOCI+IBFN
       DO JCHI=1,L1
C      -- DENSITY TERM
       DAIJ =DA(JCHI,ICHI)+DA(JCHI,ICHI)
       DBIJ =DB(JCHI,ICHI)+DB(JCHI,ICHI)
       PAIJ=PA(ICHI,JCHI,1)+PA(JCHI,ICHI,1)
       PBIJ=PA(ICHI,JCHI,IRU)+PA(JCHI,ICHI,IRU)
       VAIJ=V1A(ICHI,JCHI,1)+V1A(JCHI,ICHI,1)
       VBIJ=V1A(ICHI,JCHI,IRU)+V1A(JCHI,ICHI,IRU)
C     --ALPHA
       DUMNA=DUMA(1,1)*DAIJ+DEMA(1,1)*PAIJ+DOMA(1,1)*VAIJ
       DUMXA=DUMA(2,1)*DAIJ+DEMA(2,1)*PAIJ+DOMA(2,1)*VAIJ
       DUMYA=DUMA(3,1)*DAIJ+DEMA(3,1)*PAIJ+DOMA(3,1)*VAIJ
       DUMZA=DUMA(4,1)*DAIJ+DEMA(4,1)*PAIJ+DOMA(4,1)*VAIJ
C     --BETA
       DUMNB=DUMA(1,2)*DBIJ+DEMA(1,2)*PBIJ+DOMA(1,2)*VBIJ
       DUMXB=DUMA(2,2)*DBIJ+DEMA(2,2)*PBIJ+DOMA(2,2)*VBIJ
       DUMYB=DUMA(3,2)*DBIJ+DEMA(3,2)*PBIJ+DOMA(3,2)*VBIJ
       DUMZB=DUMA(4,2)*DBIJ+DEMA(4,2)*PBIJ+DOMA(4,2)*VBIJ
C     --TOTAL
       DUM=DUMNA+DUMNB
       DUMX=DUMXA+DUMXB
       DUMY=DUMYA+DUMYB
       DUMZ=DUMZA+DUMZB
C
C      TEMP =DUM * XYZDDC( 1,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)
C    *      +DUMX*(XYZDDC( 2,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
C    *             XYZDDC( 1,ICHI,IIPT)*XYZDDC( 2,JCHI,IIPT))
C    *      +DUMY*(XYZDDC( 3,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
C    *             XYZDDC( 1,ICHI,IIPT)*XYZDDC( 3,JCHI,IIPT))
C    *      +DUMZ*(XYZDDC( 4,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
C    *             XYZDDC( 1,ICHI,IIPT)*XYZDDC( 4,JCHI,IIPT))
       TEMPX=DUM * XYZDDC( 2,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)
     *      +DUMX*(XYZDDC( 5,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
     *             XYZDDC( 2,ICHI,IIPT)*XYZDDC( 2,JCHI,IIPT))
     *      +DUMY*(XYZDDC( 8,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
     *             XYZDDC( 2,ICHI,IIPT)*XYZDDC( 3,JCHI,IIPT))
     *      +DUMZ*(XYZDDC( 9,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
     *             XYZDDC( 2,ICHI,IIPT)*XYZDDC( 4,JCHI,IIPT))
       TEMPY=DUM * XYZDDC( 3,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)
     *      +DUMX*(XYZDDC( 8,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
     *             XYZDDC( 3,ICHI,IIPT)*XYZDDC( 2,JCHI,IIPT))
     *      +DUMY*(XYZDDC( 6,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
     *             XYZDDC( 3,ICHI,IIPT)*XYZDDC( 3,JCHI,IIPT))
     *      +DUMZ*(XYZDDC(10,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
     *             XYZDDC( 3,ICHI,IIPT)*XYZDDC( 4,JCHI,IIPT))
       TEMPZ=DUM * XYZDDC( 4,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)
     *      +DUMX*(XYZDDC( 9,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
     *             XYZDDC( 4,ICHI,IIPT)*XYZDDC( 2,JCHI,IIPT))
     *      +DUMY*(XYZDDC(10,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
     *             XYZDDC( 4,ICHI,IIPT)*XYZDDC( 3,JCHI,IIPT))
     *      +DUMZ*(XYZDDC( 7,ICHI,IIPT)*XYZDDC( 1,JCHI,IIPT)+
     *             XYZDDC( 4,ICHI,IIPT)*XYZDDC( 4,JCHI,IIPT))
C     --- USING TRANSLATIONAL INVARIANCE
       DE(1,IAT)=DE(1,IAT)-TEMPX*XYZWGT(IPT)
       DE(2,IAT)=DE(2,IAT)-TEMPY*XYZWGT(IPT)
       DE(3,IAT)=DE(3,IAT)-TEMPZ*XYZWGT(IPT)
       DE(1,KAT)=DE(1,KAT)+TEMPX*XYZWGT(IPT)
       DE(2,KAT)=DE(2,KAT)+TEMPY*XYZWGT(IPT)
       DE(3,KAT)=DE(3,KAT)+TEMPZ*XYZWGT(IPT)
       END DO
       END DO
  100 CONTINUE
C     -- DEBUG DENSITY
      RHO =RHO +RHOT *XYZWGT(IPT)
      ENDIF
      END DO
C     -- END GRID LOOP
C
      RETURN
      END
