C 11 APR 08 - MWS - REWRITE ALL GRDXX EVALUATIONS FROM **TWO TO MULTS
C  4 MAR 08 - SS  - INTERFACED EXISITING FUNCTIONALS WITH TDDFT ENERGY
C 20 AUG 07 - MC  - ADD THIRD DERIVATIVES OF FUNCTIONALS
C 29 JUN 07 - MC  - CHANGE NUMERICAL INTEGRATION ALGORHISM
C 22 DEC 06 - MC  - ADD LC EXCHANGE FUNCTIONALS
C 19 NOV 06 - MC  - NDFVWN: FIX OPEN SHELL PART OF THE FUNCTIONAL
C  6 NOV 06 - MC  - NEW MODULE FOR TDDFT'S FUNCTIONALS
C
C*MODULE TDFFUN  *DECK XCINDX
      SUBROUTINE XCINDX
      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
C
C     ---- EXCHANGE
C
C     ---- 1ST DERIVATIVES
      KRA=1
      KRB=2
      KGA=3
      KGB=4
      NXDIM(1)=4
C     ---- 2ND DERIVATIVES
      KRARA=5
      KRAGA=6
      KGAGA=7
      KRBRB=8
      KRBGB=9
      KGBGB=10
      NXDIM(2)=10
C     ---- 3RD DERIVATIVES
      KRARARA=11
      KRARAGA=12
      KRAGAGA=13
      KGAGAGA=14
      KRBRBRB=15
      KRBRBGB=16
      KRBGBGB=17
      KGBGBGB=18
      NXDIM(3)=18
C
C     ---- CORRELATION
C
C     ---- 1ST DERIVATIVES
      IRA=1
      IRB=2
      IGA=3
      IGB=4
      IGC=5
      NCDIM(1)=5
C     ---- 2ND DERIVATIVES
      IRARA=6
      IRARB=7
      IRAGA=8
      IRAGB=9
      IRAGC=10
      IRBRB=11
      IRBGA=12
      IRBGB=13
      IRBGC=14
      IGAGA=15
      IGAGB=16
      IGAGC=17
      IGBGB=18
      IGBGC=19
      IGCGC=20
      NCDIM(2)=20
C     ---- 3RD DERIVATIVES
      IRARARA=21
      IRARARB=22
      IRARAGA=23
      IRARAGB=24
      IRARAGC=25
      IRARBRB=26
      IRARBGA=27
      IRARBGB=28
      IRARBGC=29
      IRAGAGA=30
      IRAGAGB=31
      IRAGAGC=32
      IRAGBGB=33
      IRAGBGC=34
      IRAGCGC=35
      IRBRBRB=36
      IRBRBGA=37
      IRBRBGB=38
      IRBRBGC=39
      IRBGAGA=40
      IRBGAGB=41
      IRBGAGC=42
      IRBGBGB=43
      IRBGBGC=44
      IRBGCGC=45
      IGAGAGA=46
      IGAGAGB=47
      IGAGAGC=48
      IGAGBGB=49
      IGAGBGC=50
      IGAGCGC=51
      IGBGBGB=52
      IGBGBGC=53
      IGBGCGC=54
      IGCGCGC=55
      NCDIM(3)=55
      RETURN
      END
C*MODULE TDFFUN  *DECK DFXCSET
      SUBROUTINE DFXCSET(EX0,EX,EC0,EC,DRHOI_,XYZWGT,RHOMIN,EXC,ECOR,
     * ILENG,NPTGRD,IRU,NDER,IST,IEND)
C
C     CALCULATE XC-FUNCTIONAL ENERGY AND DERIVATIVES
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL LCFLAG,LRINT
C
      COMMON /CPBE0 / CPBEX
      COMMON /CX3LYP/ CPW91
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      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
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DIMENSION XYZWGT(*)
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION EX0(*)
      DIMENSION EC0(*)
      DIMENSION EX(ILENG,*)
      DIMENSION EC(ILENG,*)
      DATA ZERO/0.0D+00/
      DATA ONE/1.0D+00/
C
C     --- CLEAR
C     --EXCHANGE
      CALL VCLR(EX0,1,ILENG)
      IXDIM=NXDIM(NDER)
      CALL VCLR(EX,1,ILENG*IXDIM)
      EXC=ZERO
C     --CORRELATION
      CALL VCLR(EC0,1,ILENG)
      ICDIM=NCDIM(NDER)
      CALL VCLR(EC,1,ILENG*ICDIM)
      ECOR=ZERO
C
C     --- EXCHANGE FUNCTIONAL
C
C
C     --- SLATER EXCHANGE ---
C
      IF(NEXFG.EQ.1) THEN
         CALL NDFSLT(EXC,EX0,EX,DRHOI_,XYZWGT,CSLT,RHOMIN,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 200
      END IF
C
C     --- BECKE 88 EXCHANGE ---
C
C         BECKE IS THE ONLY ONE THAT SUPPORTS "LONG RANGE CORRECTION"
C
      IF (NEXFG.EQ.2) THEN
         IF(LCFLAG) THEN
            CALL NDFLCB88(EXC,EX0,EX,DRHOI_,XYZWGT,RHOMIN,
     *                    ILENG,NPTGRD,NDER,IRU,IST,IEND)
         ELSE
            CALL NDFB88(EXC,EX0,EX,DRHOI_,XYZWGT,CB88,RHOMIN,
     *                  ILENG,NPTGRD,NDER,IRU,IST,IEND)
         END IF
         GO TO 200
      END IF
C
C     --- PERDEW-WANG 91 EXCHANGE ---
C
      IF(NEXFG.EQ.3) THEN
         IF (NDER.EQ.3) GOTO 400
         CALL NDFPW91XF(EXC,EX0,EX,DRHOI_,XYZWGT,CPW91,RHOMIN,
     *                  ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 200
      END IF
C
C     --- PBE EXCHANGE ---
C
      IF(NEXFG.EQ.4) THEN
         IF (NDER.EQ.3) GOTO 400
         CALL NDFPBEXF(EXC,EX0,EX,DRHOI_,XYZWGT,CPBEX,RHOMIN,
     *                 ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 200
      END IF
C
C     --- GILL 96 EXCHANGE ---
C
C     IF (NEXFG.EQ.5) THEN
C        GO TO 200
C     END IF
C
C     --- PARAMETER FREE EXCHANGE ---
C
C     IF (NEXFG.EQ.6) THEN
C        GO TO 200
C     END IF
C
C     --- OPTX EXCHANGE ---
C
C     IF (NEXFG.EQ.7) THEN
C        GO TO 200
C     END IF
C
C     --- B3 HYBRID'S FUNCTIONAL PARTS
C
      IF (NEXFG.EQ.-1) THEN
         CALL NDFSLT(EXC,EX0,EX,DRHOI_,XYZWGT,CSLT,RHOMIN,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
         CALL NDFB88(EXC,EX0,EX,DRHOI_,XYZWGT,CB88,RHOMIN,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 200
      END IF
C
C     --- X3LYP EXCHANGE ---
C
      IF (NEXFG.EQ.-2) THEN
         IF (NDER.EQ.3) GOTO 400
         CALL NDFSLT(EXC,EX0,EX,DRHOI_,XYZWGT,CSLT,RHOMIN,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
         CALL NDFB88(EXC,EX0,EX,DRHOI_,XYZWGT,CB88,RHOMIN,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
         CALL NDFPW91XF(EXC,EX0,EX,DRHOI_,XYZWGT,CPW91,RHOMIN,
     *                  ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 200
      END IF
C
      WRITE(IW,*) 'TDDFUN: EXCHANGE FUNCTIONAL IS NOT SUPPORTED'
      CALL ABRT
      STOP
C
C     -- CORRELATION FUNCTIONAL
C
  200 CONTINUE
C
C     -- LOOSENING THREHOLD FOR THIRD DERIVATIVE OF LYP FUNCTIONAL
      RHOMINL=RHOMIN
      IF(NDER.EQ.3) RHOMINL=1.0D-10
C
      IF(NCORFG.EQ.0) THEN
C        WRITE(IW,*) 'NO CORRELATION FUNCTIONAL USED'
         GO TO 300
      END IF
C
C     --- VWN5 CORRELATION FUNCTIONAL ---
C
      IF(NCORFG.EQ.1) THEN
         CALL NDFVWN(ECOR,EC0,EC,DRHOI_,XYZWGT,CVWN,RHOMIN,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
C     --- LEE-YANG-PARR CORRELATION FUNCTIONAL ---
C
      IF(NCORFG.EQ.2) THEN
         CALL NDFLYP(ECOR,EC0,EC,DRHOI_,XYZWGT,CLYP,RHOMINL,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
C     --- PERDEW-WANG 91 CORRELATION FUNCTIONAL ---
C
      IF(NCORFG.EQ.3) THEN
         IF (NDER.EQ.3) GOTO 400
         CALL NDFPW91CF(ECOR,EC0,EC,DRHOI_,XYZWGT,ONE,RHOMIN,
     *                  ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
C     --- PBE CORRELATION FUNCTIONAL ---
C
      IF(NCORFG.EQ.4) THEN
         IF (NDER.EQ.3) GOTO 400
         CALL NDFPBECF(ECOR,EC0,EC,DRHOI_,XYZWGT,ONE,RHOMIN,
     *                 ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
C     --- VWN FORMULA 1 CORRELATION FUNCTIONAL ---
C
      IF(NCORFG.EQ.5) THEN
         IF (NDER.EQ.3) GOTO 400
         CALL NDFVWN1(ECOR,EC0,EC,DRHOI_,XYZWGT,CVWN,RHOMIN,
     *                ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
C     --- PW91 LDA CORRELATION FUNCTIONAL ---
C
      IF(NCORFG.EQ.8) THEN
         IF (NDER.EQ.3) GOTO 400
         CALL NDFPW91LDA(ECOR,EC0,EC,DRHOI_,XYZWGT,ONE,RHOMIN,
     *                   ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
C     --- PERDEW 86 CORRELATION FUNCTIONAL ---
C
      IF(NCORFG.EQ.9) THEN
         IF (NDER.EQ.3) GOTO 400
         CALL NDFP86CF(ECOR,EC0,EC,DRHOI_,XYZWGT,ONE,RHOMIN,
     *                 ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
C     --- PERDEW-ZUNGER 81 CORRELATION FUNCTIONAL ---
C
      IF(NCORFG.EQ.10) THEN
         IF (NDER.EQ.3) GOTO 400
         CALL NDFPZ81CF(ECOR,EC0,EC,DRHOI_,XYZWGT,ONE,RHOMIN,
     *                  ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
C     --- B3LYP = B3LYP5 CORRELATION FUNCTIONAL ---
C
      IF(NCORFG.EQ.-1) THEN
C            -- LYP/VWN5 PART OF B3LYP
         CALL NDFVWN(ECOR,EC0,EC,DRHOI_,XYZWGT,CVWN,RHOMIN,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
         CALL NDFLYP(ECOR,EC0,EC,DRHOI_,XYZWGT,CLYP,RHOMINL,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
C     --- B3LYP1 CORRELATION FUNCTIONAL ---
C
      IF(NCORFG.EQ.-2) THEN
C            -- LYP/VWN1 PART OF B3LYP
         CALL NDFVWN1(ECOR,EC0,EC,DRHOI_,XYZWGT,CVWN,RHOMIN,
     *                ILENG,NPTGRD,NDER,IRU,IST,IEND)
         CALL NDFLYP(ECOR,EC0,EC,DRHOI_,XYZWGT,CLYP,RHOMINL,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
C     --- PERDEW-WANG 91 CORRELATION FUNCTIONAL ---
C     --- UTILIZING NWCHEM CONSTANTS FOR PW91LDA ---
C
      IF(NCORFG.EQ.-3) THEN
         IF (NDER.EQ.3) GOTO 400
         CALL NDFNWPW91CF(ECOR,EC0,EC,DRHOI_,XYZWGT,ONE,RHOMIN,
     *                    ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
C     --- ONE-PARAMETER PROGRESSIVE CORRELATION FUNCTIONAL ---
C
C     THIS FUNCTIONAL WAS SWITCHED FROM NCORFG.EQ.4 TO NCORFG.EQ.7
C     IN ORDER TO FOLLOW THE NUMBER SCHEME FROM ROUTINE CALCEXC IN
C     DFTEXC.SRC
C
      IF(NCORFG.EQ.7) THEN
         CALL NDFOP(ECOR,EC0,EC,DRHOI_,XYZWGT,RHOMIN,
     *              ILENG,NPTGRD,NDER,IRU,IST,IEND)
         GO TO 300
      END IF
C
      WRITE(IW,*) 'TDDFUN: CORRELATION FUNCTIONAL IS NOT SUPPORTED'
      CALL ABRT
      STOP
C
  400 CONTINUE
      WRITE(IW,*) '3RD FUNCTIONAL DERIVATIVES ARE NOT YET AVAILABLE'
      CALL ABRT
      STOP
C
  300 CONTINUE
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFSLT
      SUBROUTINE NDFSLT(EXC,EX0,EX,DRHOI_,XYZWGT,SCALE,RHOMIN,
     * ILENG,NPTGRD,NDER,IRU,IST,IEND)
      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
      DIMENSION EX0(*),EX(ILENG,*),DRHOI_(NPTGRD,4,*)
      DIMENSION RHOS(2)
      DIMENSION XYZWGT(*)
      DATA ZERO      /0.0D+00/
      DATA ONE       /1.0D+00/
      DATA TWO       /2.0D+00/
      DATA THREE     /3.0D+00/
      DATA FOUR      /4.0D+00/
      DATA FIVE      /5.0D+00/
C
      ESLT=ZERO
      RHO=ZERO
      DO IPT=IST,IEND
      IIPT=IPT-IST+1
       RHOS(KRA)=DRHOI_(IPT,1,1)
      IF(IRU.EQ.1) THEN
       RHOS(KRB)=RHOS(KRA)
      ELSE
       RHOS(KRB)=DRHOI_(IPT,1,IRU)
      ENDIF
C     -- SPIN LOOP
      DO I=KRA,KRB
      IF(I.EQ.1) THEN
       IR=KRA
       IRR=KRARA
       IRRR=KRARARA
      ELSE
       IR=KRB
       IRR=KRBRB
       IRRR=KRBRBRB
      ENDIF
      RHOI=RHOS(I)
      IF(RHOI.GT.RHOMIN) THEN
C
C     ----- SLATER POTENTIAL
C
      PI  =FOUR*ATAN(ONE)
      FACE=-(THREE/TWO)*(THREE/PI/FOUR)**(ONE/THREE)
      FACP=FOUR/THREE*FACE
      FACP2=ONE/THREE*FACP
      FACP3=-TWO/THREE*FACP2
C
      E0=SCALE*FACE*RHOI**(FOUR/THREE)
      EX0(IIPT)=EX0(IIPT)+E0
C    ---- 1ST DERIVATIVE
      EX(IIPT,IR)=EX(IIPT,IR)+SCALE*FACP*RHOI**(ONE/THREE)
C     -- ENERGY & DENSITY CHECK
      ESLT=ESLT+E0*XYZWGT(IPT)
      RHO   =RHO   +RHOI *XYZWGT(IPT)
C
      IF (NDER.EQ.1) GOTO 100
C    ---- 2ND DERIVATIVE
      EX(IIPT,IRR)=EX(IIPT,IRR)+SCALE*FACP2*RHOI**(-TWO/THREE)
      IF (NDER.EQ.2) GOTO 100
C     ---- 3RD DERIVATIVE
      EX(IIPT,IRRR)=EX(IIPT,IRRR)+SCALE*FACP3*RHOI**(-FIVE/THREE)
      ENDIF
  100 CONTINUE
      END DO
C     -- END SPIN LOOP
      END DO
C     -- END GRID LOOP
      EXC=EXC+ESLT
C
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFB88
      SUBROUTINE NDFB88(EXC,EX0,EX,DRHOI_,XYZWGT,SCALE,RHOMIN,
     * ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION K0,KX0
      LOGICAL DBUG
      COMMON /XDERIX/ NXDIM(3),KRA,KRB,KGA,KGB,KRARA,KRAGA,KGAGA,
     *                KRBRB,KRBGB,KGBGB,KRARARA,KRARAGA,KRAGAGA,KGAGAGA,
     *                KRBRBRB,KRBRBGB,KRBGBGB,KGBGBGB
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION EX0(*)
      DIMENSION EX(ILENG,*)
      DIMENSION XYZWGT(*)
      DIMENSION RHOS(2)
      DIMENSION DRHOS(2)
      DIMENSION RHOCHECK(2)
      DIMENSION EXCS(2)
      DATA ZERO/0.0D+00/
      DATA HALF/0.5D+00/
      DATA QUATER/0.25D+00/
      DATA ONEITH/0.125D+00/
      DATA ONE/1.0D+00/
      DATA TWO/2.0D+00/
      DATA THREE/3.0D+00/
      DATA FOUR/4.0D+00/
      DATA SIX/6.0D+00/
      DATA EIGHT/8.0D+00/
      DATA SEVEN/7.0D+00/
      DATA XNINE/9.0D+00/
      DATA FORTN/14.0D+00/
      DBUG=.FALSE.
      ONETHD=ONE/THREE
      FORTHD=FOUR/THREE
      EGHT27=EIGHT/27.0D+00
      TWONTH=TWO/XNINE
C
      CB88=0.0042D+00
      PI  =ATAN(ONE)*FOUR
      K0=-1.5D+00*(0.75D+00/PI)**ONETHD
      IF(DBUG) THEN
      ENDIF
C     -- CLEAR
      DO I=1,2
      RHOCHECK(I)=ZERO
      EXCS(I)=ZERO
      END DO
C     -- GRID LOOP
      DO IPT=IST,IEND
      IIPT=IPT-IST+1
      RHOS(KRA)=DRHOI_(IPT,1,1)
      DRHOS(KRA)=DRHOI_(IPT,2,1)**2+DRHOI_(IPT,3,1)**2
     * +DRHOI_(IPT,4,1)**2
      IF(IRU.EQ.1) THEN
       RHOS(KRB)=RHOS(KRA)
       DRHOS(KRB)=DRHOS(KRA)
      ELSE
       RHOS(KRB)=DRHOI_(IPT,1,IRU)
       DRHOS(KRB)=DRHOI_(IPT,2,IRU)**2+DRHOI_(IPT,3,IRU)**2
     * +DRHOI_(IPT,4,IRU)**2
      ENDIF
C     -- SPIN LOOP
      DO I=KRA,KRB
      IF(I.EQ.1) THEN
      IR=KRA
      IG=KGA
      IRR=KRARA
      IRG=KRAGA
      IGG=KGAGA
      IRRR=KRARARA
      IRRG=KRARAGA
      IRGG=KRAGAGA
      IGGG=KGAGAGA
      ELSE
      IR=KRB
      IG=KGB
      IRR=KRBRB
      IRG=KRBGB
      IGG=KGBGB
      IRRR=KRBRBRB
      IRRG=KRBRBGB
      IRGG=KRBGBGB
      IGGG=KGBGBGB
      ENDIF
      RHO=RHOS(I)
      GM=DRHOS(I)
      IF (RHO.GT.RHOMIN) THEN
C     -- FOR BECKE88 DERIVATIVES
      RHO13=RHO**(ONE/THREE)
      RHO43=RHO13 * RHO
      GM12  =SQRT(GM)
      X     =GM12/RHO43
      X2     = X*X
      X2P1   = ONE+X2
      X2P1_12= SQRT(X2P1)
      ASINH  = LOG(X+X2P1_12)
      AA  = CB88 * X2
      BB  = ONE + SIX * CB88*X * ASINH
      FX     = - AA/BB
      C6  = CB88*SIX
      DAA = TWO*CB88*X
      DBB = C6 * (ASINH + X/X2P1_12)
      DD  = DAA*BB-AA*DBB
      BB2 = BB*BB
      DFDX  = -DD/BB2
      KX0=-TWO*(K0 + FX)
C     -- B88 ENERGY
      EE0=-SCALE*HALF* KX0 * RHO43
      EX0(IIPT)=EX0(IIPT)-SCALE*HALF* KX0 * RHO43
C     -----------------------
      TMP2=K0+FX-DFDX*X
C     -- B88 1ST DERIVATIVES
      EX(IIPT,IR)=EX(IIPT,IR)+SCALE*FORTHD * TMP2 * RHO13
      EX(IIPT,IG)=EX(IIPT,IG)+SCALE* HALF * DFDX /GM12
C     -- ENERGY & DENSITY CHECK
      EXCS(I)=EXCS(I)+EE0*XYZWGT(IPT)
      RHOCHECK(I)=RHOCHECK(I)+RHO*XYZWGT(IPT)
      IF (NDER.EQ.1) GOTO 100
C     -----------------------------
      X2P1_32= X2P1_12 * X2P1
      D2AA=TWO*CB88
      D2BB=C6*(TWO/X2P1_12 -X2/X2P1_32)
      DD1 =D2AA*BB-D2BB*AA
      EE  =DD1*BB-TWO*DD*DBB
      BB3 =BB2*BB
      D2FDX=-EE/BB3
      X2    =X*X
      RHO23M=RHO13/RHO
      RHO73 =RHO43*RHO
      GM32  =GM12 *GM
      TMP2=TMP2 + FOUR*D2FDX*X2
      TMP =D2FDX*X-DFDX
C     -- B88 2ND DERIVATIVES
      EX(IIPT,IRR)=EX(IIPT,IRR)+SCALE*FOUR/XNINE * TMP2 * RHO23M
      EX(IIPT,IRG)=EX(IIPT,IRG)-SCALE*TWO/THREE * D2FDX / RHO73
      EX(IIPT,IGG)=EX(IIPT,IGG)+SCALE* QUATER * TMP / GM32
      IF (NDER.EQ.2) GOTO 100
C     ---------------------------------
      X2P1_52= X2P1_32 * X2P1
      D3BB =C6*(-FOUR*X/X2P1_32 +THREE*X2*X/X2P1_52)
      DD2 = D2AA*DBB - DAA*D2BB - AA*D3BB
      DEE = DD2*BB -DD1*DBB -TWO*DD*D2BB
      FF  = DEE*BB -THREE*EE*DBB
      BB4 = BB3*BB
      D3FDX  = -FF/BB4
      X3    = X2*X
      RHO53M= RHO23M/RHO
      RHO103= RHO73*RHO
      GM52  = GM32 * GM
      TMP   = D3FDX*X2 -TMP * THREE
      TMP2  = TMP2 + EIGHT*D3FDX*X3 + FORTN*D2FDX*X2
C     -- B88 3RD DERIVATIVES
      EX(IIPT,IRRR)=EX(IIPT,IRRR)-SCALE*EGHT27*TMP2*RHO53M
      EX(IIPT,IRRG)=EX(IIPT,IRRG)+SCALE*TWONTH*(FOUR*D3FDX*X
     * +SEVEN*D2FDX)/RHO103
      EX(IIPT,IRGG)=EX(IIPT,IRGG)-SCALE*ONETHD* D3FDX*X /RHO73/GM
      EX(IIPT,IGGG)=EX(IIPT,IGGG)+SCALE*ONEITH * TMP / GM52
C     -------------------------------------------------
      ENDIF
  100 CONTINUE
      END DO
C     -- SPIN LOOP END
      END DO
C     -- GRID LOOP END
      RHOT=RHOCHECK(1)+RHOCHECK(2)
      EB88=EXCS(1)+EXCS(2)
      EXC=EXC+EB88
C
      IF (DBUG) THEN
         WRITE(6,*)'B88 EXCHANGE ENERGY IN DERB8D=',EB88
         WRITE(6,*)'INTEGRATED DENSITY IN DERB88D=',RHOT
      ENDIF
C
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFLCB88
      SUBROUTINE NDFLCB88(EXC,EX0,EX,DRHOI_,XYZWGT,RHOMIN,
     * ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION K0,KX0,KX12
      LOGICAL DBUG,LCFLAG,LRINT
      PARAMETER (TAYTHR=1.0D-2)
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      COMMON /DFTEXC/ PI,CSLT,CB88H,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /XDERIX/ NXDIM(3),KRA,KRB,KGA,KGB,KRARA,KRAGA,KGAGA,
     *                KRBRB,KRBGB,KGBGB,KRARARA,KRARAGA,KRAGAGA,KGAGAGA,
     *                KRBRBRB,KRBRBGB,KRBGBGB,KGBGBGB
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION EX0(*)
      DIMENSION EX(ILENG,*)
      DIMENSION XYZWGT(*)
      DIMENSION RHOS(2)
      DIMENSION DRHOS(2)
      DIMENSION RHOCHECK(2)
      DIMENSION EXCS(2)
      DIMENSION EB88(2)
      DATA ZERO/0.0D+00/
      DATA HALF/0.5D+00/
      DATA QUATER/0.25D+00/
      DATA ONEITH/0.125D+00/
      DATA ONE/1.0D+00/
      DATA TWO/2.0D+00/
      DATA THREE/3.0D+00/
      DATA FOUR/4.0D+00/
      DATA SIX/6.0D+00/
      DATA SEVEN/7.0D+00/
      DATA EIGHT/8.0D+00/
      DATA XNINE/9.0D+00/
      DATA FORTN/14.0D+00/
      DBUG=.FALSE.
      ONETHD=ONE/THREE
      FORTHD=FOUR/THREE
      EGHT27=EIGHT/27.0D+00
      TWONTH=TWO/XNINE
C
C     --DEBUG
C     EMU=ZERO
      CB88=0.0042D+00
      PI  =ATAN(ONE)*FOUR
      PI12=SQRT(PI)
      CA_S=EMU/(SIX*PI12)
      K0=-1.5D+00*(0.75D+00/PI)**ONETHD
C     -- CLEAR
      DO I=1,2
      RHOCHECK(I)=ZERO
      EXCS(I)=ZERO
      EB88(I)=ZERO
      END DO
C     -- GRID LOOP
      DO IPT=IST,IEND
      IIPT=IPT-IST+1
      RHOS(KRA)=DRHOI_(IPT,1,1)
      DRHOS(KRA)=DRHOI_(IPT,2,1)**2+DRHOI_(IPT,3,1)**2
     * +DRHOI_(IPT,4,1)**2
      IF(IRU.EQ.1) THEN
       RHOS(KRB)=RHOS(KRA)
       DRHOS(KRB)=DRHOS(KRA)
      ELSE
       RHOS(KRB)=DRHOI_(IPT,1,IRU)
       DRHOS(KRB)=DRHOI_(IPT,2,IRU)**2+DRHOI_(IPT,3,IRU)**2
     * +DRHOI_(IPT,4,IRU)**2
      ENDIF
C     -- SPIN LOOP
      DO I=KRA,KRB
      IF(I.EQ.1) THEN
      IR=KRA
      IG=KGA
      IRR=KRARA
      IRG=KRAGA
      IGG=KGAGA
      IRRR=KRARARA
      IRRG=KRARAGA
      IRGG=KRAGAGA
      IGGG=KGAGAGA
      ELSE
      IR=KRB
      IG=KGB
      IRR=KRBRB
      IRG=KRBGB
      IGG=KGBGB
      IRRR=KRBRBRB
      IRRG=KRBRBGB
      IRGG=KRBGBGB
      IGGG=KGBGBGB
      ENDIF
      RHO=RHOS(I)
      GM=DRHOS(I)
      IF (RHO.GT.RHOMIN) THEN
C     -- FOR BECKE88 DERIVATIVES
      RHO13=RHO**(ONE/THREE)
      RHO43=RHO13 * RHO
      GM12  =SQRT(GM)
      X     =GM12/RHO43
      X2     = X*X
      X2P1   = ONE+X2
      X2P1_12= SQRT(X2P1)
      ASINH  = LOG(X+X2P1_12)
      AA  = CB88 * X2
      BB  = ONE + SIX * CB88*X * ASINH
      FX     = - AA/BB
      C6  = CB88*SIX
      DAA = TWO*CB88*X
      DBB = C6 * (ASINH + X/X2P1_12)
      DD  = DAA*BB-AA*DBB
      BB2 = BB*BB
      DFDX  = -DD/BB2
      KX0=-TWO*(K0 + FX)
C     -- B88 ENERGY
      E0= -HALF* KX0 * RHO43
C     -----------------------
      TMP2=K0+FX-DFDX*X
C     -- B88 1ST DERIVATIVES
      DEDR=  FORTHD * TMP2 * RHO13
      DEDG=  HALF * DFDX /GM12
C     -----------------------------
      X2P1_32= X2P1_12 * X2P1
      D2AA=TWO*CB88
      D2BB=C6*(TWO/X2P1_12 -X2/X2P1_32)
      DD1 =D2AA*BB-D2BB*AA
      EE  =DD1*BB-TWO*DD*DBB
      BB3 =BB2*BB
      D2FDX=-EE/BB3
      X2    =X*X
      RHO23M=RHO13/RHO
      RHO73 =RHO43*RHO
      GM32  =GM12 *GM
      TMP2=TMP2 + FOUR*D2FDX*X2
      TMP =D2FDX*X-DFDX
C     -- B88 2ND DERIVATIVES
      DEDRR= FOUR/XNINE * TMP2 * RHO23M
      DEDRG=-TWO/THREE * D2FDX / RHO73
      DEDGG= QUATER * TMP / GM32
C     ---------------------------------
      X2P1_52= X2P1_32 * X2P1
      D3BB =C6*(-FOUR*X/X2P1_32 +THREE*X2*X/X2P1_52)
      DD2 = D2AA*DBB - DAA*D2BB - AA*D3BB
      DEE = DD2*BB -DD1*DBB -TWO*DD*D2BB
      FF  = DEE*BB -THREE*EE*DBB
      BB4 = BB3*BB
      D3FDX  = -FF/BB4
      X3    = X2*X
      RHO53M= RHO23M/RHO
      RHO103= RHO73*RHO
      GM52  = GM32 * GM
      TMP   = D3FDX*X2 -TMP * THREE
      TMP2  = TMP2 + EIGHT*D3FDX*X3 + FORTN*D2FDX*X2
C     -- B88 3RD DERIVATIVES
      DEDRRR=-EGHT27*TMP2*RHO53M
      DEDRRG= TWONTH*(FOUR*D3FDX*X+SEVEN*D2FDX)/RHO103
      DEDRGG=-ONETHD* D3FDX*X /RHO73/GM
      DEDGGG=ONEITH * TMP / GM52
C     -------------------------------------------------
C     -- FOR LONG-RANGE CORRECTION TERM
      KX12=SQRT(KX0)
      A_S = CA_S*KX12/RHO13
      A_S2= A_S**2
C     -- ERROR FUNCTION
      XERF = DFTERF(HALF/A_S)
      TMP=-QUATER/A_S2
      IF (ABS(TMP).LT.TAYTHR) THEN
        B_S = TAYEXP(TMP)
      ELSE
        B_S = EXP(TMP)-ONE
      ENDIF
      C_S = TWO * A_S2 * B_S + HALF
      D_S = PI12*XERF + TWO*A_S*(B_S-C_S)
      R_S =EIGHT/THREE * A_S * D_S
      F0 = ONE-R_S
C     ** ENERGY **
      EX0(IIPT)=EX0(IIPT)+E0 * F0
C     *************
C     -- G-TERM
      G=SQRT(-E0)
      R2=SQRT(TWO)
C     -- FOR FIRST DERIVATIVES
      DGR=-HALF/G*DEDR
      DGG=-HALF/G*DEDG
      DAR=ONE/RHO*(CA_S*R2*DGR-A_S)
      DAG=ONE/RHO*(CA_S*R2*DGG)
      DRDA=EIGHT/THREE*(D_S-SIX*A_S*C_S)
      DFR=-DRDA*DAR
      DFG=-DRDA*DAG
C     ** FIRST DERIVATIVES **/
      EX(IIPT,IR)=E0*DFR+DEDR*F0
      EX(IIPT,IG)=E0*DFG+DEDG*F0
C     -- ENERGY & DENSITYCHECK
      EXCS(I)=EXCS(I)+E0 * F0*XYZWGT(IPT)
      EB88(I)=EB88(I)+E0*XYZWGT(IPT)
      RHOCHECK(I)=RHOCHECK(I)+RHO*XYZWGT(IPT)
      IF (NDER.EQ.1) GOTO 100
C     **************************
C     -- FOR SECOND DERIVATIVES
      DGRR=-ONE/G*(DGR*DGR+HALF*DEDRR)
      DGRG=-ONE/G*(DGR*DGG+HALF*DEDRG)
      DGGG=-ONE/G*(DGG*DGG+HALF*DEDGG)
      DARR=ONE/RHO*(CA_S*R2*DGRR-DAR-DAR)
      DARG=ONE/RHO*(CA_S*R2*DGRG-DAG)
      DAGG=ONE/RHO*(CA_S*R2*DGGG)
      DRDAA=-64.0D+00*C_S-16.0D+00*B_S
      DFRR=-DRDAA*DAR*DAR-DRDA*DARR
      DFRG=-DRDAA*DAR*DAG-DRDA*DARG
      DFGG=-DRDAA*DAG*DAG-DRDA*DAGG
C     ** SECOND DERIVATIVES **/
      EX(IIPT,IRR)=E0*DFRR+TWO*DEDR*DFR+DEDRR*F0
      EX(IIPT,IRG)=E0*DFRG+DEDR*DFG+DEDG*DFR+DEDRG*F0
      EX(IIPT,IGG)=E0*DFGG+TWO*DEDG*DFG+DEDGG*F0
      IF (NDER.EQ.2) GOTO 100
C     ***********************************************
C     -- FOR THIRD DERIVATIVES
      DGRRR=-ONE/G*(DGRR*DGR+DGRR*DGR+DGRR*DGR+HALF*DEDRRR)
      DGRRG=-ONE/G*(DGRG*DGR+DGRG*DGR+DGRR*DGG+HALF*DEDRRG)
      DGRGG=-ONE/G*(DGGG*DGR+DGRG*DGG+DGRG*DGG+HALF*DEDRGG)
      DGGGG=-ONE/G*(DGGG*DGG+DGGG*DGG+DGGG*DGG+HALF*DEDGGG)
      DARRR=ONE/RHO*(CA_S*R2*DGRRR-DARR-DARR-DARR)
      DARRG=ONE/RHO*(CA_S*R2*DGRRG-DARG-DARG)
      DARGG=ONE/RHO*(CA_S*R2*DGRGG-DAGG)
      DAGGG=ONE/RHO*(CA_S*R2*DGGGG)
C     DUM=(B_S+ONE)/A_S2
C     DRDAAA=-EIGHT/A_S*(16.0D+00*C_S+EIGHT*B_S+DUM)
      DCDA=FOUR*A_S*B_S+ONE/A_S*(B_S+ONE)
      DBDA=HALF/(A_S**3)*(B_S+ONE)
      DRDAAA=-64.0D+00*DCDA-16.0D+00*DBDA
      DUMRRR=DARR*DAR+DARR*DAR+DARR*DAR
      DFRRR=-DRDAAA*DAR*DAR*DAR-DRDAA*DUMRRR-DRDA*DARRR
      DUMRRG=DARG*DAR+DARG*DAR+DARR*DAG
      DFRRG=-DRDAAA*DAR*DAR*DAG-DRDAA*DUMRRG-DRDA*DARRG
      DUMRGG=DAGG*DAR+DARG*DAG+DARG*DAG
      DFRGG=-DRDAAA*DAR*DAG*DAG-DRDAA*DUMRGG-DRDA*DARGG
      DUMGGG=DAGG*DAG+DAGG*DAG+DAGG*DAG
      DFGGG=-DRDAAA*DAG*DAG*DAG-DRDAA*DUMGGG-DRDA*DAGGG
C     ** THIRD DERIVATIVES **/
      EX(IIPT,IRRR)=DEDRRR*F0+THREE*(DEDRR*DFR+DEDR*DFRR)+E0*DFRRR
      EX(IIPT,IRRG)=DEDRRG*F0+DEDRR*DFG+DEDRG*DFR+DEDRG*DFR+
     * DEDR*DFRG+DEDR*DFRG+DEDG*DFRR+E0*DFRRG
      EX(IIPT,IRGG)=DEDRGG*F0+DEDRG*DFG+DEDGG*DFR+DEDRG*DFG+
     * DEDR*DFGG+DEDG*DFRG+DEDG*DFRG+E0*DFRGG
      EX(IIPT,IGGG)=DEDGGG*F0+THREE*(DEDGG*DFG+DEDG*DFGG)+E0*DFGGG
C ******************************************************************
      ENDIF
  100 CONTINUE
      END DO
C     -- SPIN LOOP END
      END DO
C     -- GRID LOOP END
      RHOT=RHOCHECK(1)+RHOCHECK(2)
      ELCB88=EXCS(1)+EXCS(2)
      EB88_NONLC=EB88(1)+EB88(2)
      EXC=EXC+ELCB88
C
      IF (DBUG) THEN
         WRITE(6,*)'LC-B88 EXCHANGE ENERGY IN DERLCB88=',ELCB88
         WRITE(6,*)'ALPHA SPIN EXCHANGE ENERGY=',EXCS(1)
         WRITE(6,*)'BETA SPIN EXCHANGE ENERGY=',EXCS(2)
         WRITE(6,*)'NON_LC-B88 EXCHANGE ENERGY=',EB88_NONLC
         WRITE(6,*)'INTEGRATED DENSITY IN DERLCB88=',RHOT
         WRITE(6,*)'ALPHA DENSITY=',RHOCHECK(1)
         WRITE(6,*)'BETA DENSITY=',RHOCHECK(2)
         WRITE(6,*)'IRU AND NDER=',IRU,NDER
      ENDIF
C
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFOP
      SUBROUTINE NDFOP(ECOR,EC0,EC,DRHOI_,XYZWGT,RHOMIN,
     * ILENG,NPTGRD,NDER,IRU,IST,IEND)
C
C     CODE DESCRIPTION BY :  MUNEAKI KAMIYA 2004/09/12
C
C      ONE-PARAMETER PROGRESSIVE CORRELATION FUNCTIONAL
C
C                T. TSUNEDA AND K. HIRAO
C         J. CHEM. PHYS. 110, 10664-10678 (1999)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION KX0(2)
      DOUBLE PRECISION KX(3,2)
      PARAMETER(C1=1.5214D+00,C2=0.5764D+00,C3=1.1284D+00,C4=0.3183D+00)
      LOGICAL DBUG
      COMMON /DFTEXC/ PI,CSLT,CCB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      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 DRHOI_(NPTGRD,4,IRU)
      DIMENSION EC0(*)
      DIMENSION EC(ILENG,*)
      DIMENSION XYZWGT(*)
      DIMENSION RHO(2),RHO13(2),RHO43(2),GM(2),GM12(2),
     * X(2)
C     ---- RHO13KX
      DIMENSION RK(2),FF(2)
      DIMENSION RHO23(2),RHO53(2),RHO83(2),GM2(2),GM3(2)
      DIMENSION DFF(2) ,DRK_R(2)   ,DRK_G(2)
      DIMENSION D2FF(2),D2RK_RR(2) ,D2RK_RG(2) ,D2RK_GG(2)
      DIMENSION D3FF(2),D3RK_RRR(2),D3RK_RRG(2),D3RK_RGG(2),D3RK_GGG(2)
C     ---- BETA
      DIMENSION TMP(55)
C
      DATA ZERO      /0.0D+00/
      DATA PT5       /0.5D+00/
      DATA ONE       /1.0D+00/
      DATA TWO       /2.0D+00/
      DATA THREE     /3.0D+00/
      DATA FOUR      /4.0D+00/
      DATA FIVE      /5.0D+00/
      DATA SIX       /6.0D+00/
      DATA EIGHT     /8.0D+00/
      DATA XNINE     /9.0D+00/
      DBUG=.FALSE.
C
C     --- COEFFICIENTS
C
C     --1ST
      CK_R=ONE/THREE
      CK_G=PT5
C     --2ND
      CK_RR=-TWO/(THREE*THREE)
      CK_RG=ONE/SIX
      CK_GG=PT5*PT5
C     --3RD
      CK_RRR= TWO/(THREE**3)
      CK_RRG=-ONE/(THREE*THREE)
      CK_RGG=-ONE/(THREE*FOUR)
      CK_GGG= PT5*PT5
C
C     ---- GRID LOOP
C
      RHOCHECK=ZERO
      EOP=ZERO
      DO IPT=IST,IEND
      IIPT=IPT-IST+1
      RHO(IRA)=DRHOI_(IPT,1,1)
      GM(IRA) =DRHOI_(IPT,2,1)**2+DRHOI_(IPT,3,1)**2
     * +DRHOI_(IPT,4,1)**2
       IF(IRU.EQ.1) THEN
       RHO(IRB)=RHO(IRA)
       GM(IRB)=GM(IRA)
       ELSE
       RHO(IRB)=DRHOI_(IPT,1,IRU)
       GM(IRB) =DRHOI_(IPT,2,IRU)**2+DRHOI_(IPT,3,IRU)**2
     * +DRHOI_(IPT,4,IRU)**2
       ENDIF
       RHOTOT=RHO(IRA)+RHO(IRB)
C
       IF(RHOTOT.GT.RHOMIN*TWO) THEN
C
       DO I=1,2
       RHO13(I)=RHO(I)**CK_R
       RHO43(I)=RHO13(I) * RHO(I)
       GM12(I)  =SQRT(GM(I))
       X(I) =GM12(I)/RHO43(I)
       END DO
C
C     ---- GENERATE KX0(2),KX(3,2) FROM X(2)
C     --CLEAR KX
      DO I=1,2
      KX0(I)=ZERO
      END DO
      DO I=1,2
      DO J=1,3
      KX(J,I)=ZERO
      END DO
      END DO
C
      PI=ATAN(ONE)*FOUR
      FACE=-THREE/TWO*((THREE/FOUR)/PI)**(ONE/THREE)
C
      IF (NEXFG.EQ.1) THEN
C     -- FOR LDA
       QOP    =2.5654D0
       DO I=1,2
        KX0(I)=FACE
       END DO
      ELSEIF (NEXFG.EQ.2) THEN
C     -- FOR BECKE88
       CB88=0.0042D+00
       QOP =2.3670D+00
       DO I=1,2
        X2   = X(I)*X(I)
        X2P1   = ONE+X2
        X2P1_12= SQRT(X2P1)
        ASINH  = LOG(X(I)+X2P1_12)
C     --KX0(I):OK
        AA  = CB88 * X2
        BB  = ONE + SIX * CB88*X(I)* ASINH
        KX0(I)     = FACE - AA/BB
C     --KX(1,I)
        C6  = CB88*SIX
        DAA = TWO*CB88*X(I)
        DBB = C6 * (ASINH + X(I)/X2P1_12)
        DD = DAA*BB-AA*DBB
        BB2 = BB*BB
        KX(1,I)  = -DD/BB2
C     --KX(2,I)
        X2P1_32= X2P1_12 * X2P1
        D2AA=TWO*CB88
        D2BB=C6*(TWO/X2P1_12 -X2/X2P1_32)
        DD1 =D2AA*BB-D2BB*AA
        EE  =DD1*BB-2.0D0*DD*DBB
        BB3 =BB2*BB
        KX(2,I)=-EE/BB3
C     --KX(3,I)
        X2P1_52= X2P1_32 * X2P1
        D3BB =C6*(-FOUR*X(I)/X2P1_32 +THREE*X2*X(I)/X2P1_52)
        DD2 = D2AA*DBB - DAA*D2BB - AA*D3BB
        DEE = DD2*BB -DD1*DBB -TWO*DD*D2BB
        FFK  = DEE*BB -THREE*EE*DBB
        BB4 = BB3*BB
        KX(3,I)  = -FFK/BB4
       END DO
      ELSEIF (NEXFG.EQ.5) THEN
C     -- FOR G96
       GILL=1/137.0D+00
       QOP=2.3638D+00
       DO I=1,2
        X12=DSQRT(X(I))
        X32=X(I)*X12
C     --KX0(I),KX(1-3,I)
        KX0(I)=-X32*GILL+FACE
        KX(1,I)=-X12*1.5D+00*GILL
        KX(2,I)=-0.75D+00*GILL/X12
        KX(3,I)=0.375D+00*GILL/X32
       END DO
      ENDIF
C     --KX*(-TWO)
      DO I=1,2
      KX0(I)=-KX0(I)*TWO
       DO J=1,3
       KX(J,I)=-KX(J,I)*TWO
       END DO
      END DO
C     ---- GEN KX0 & KX END ----
C
C     ---- CORRELATION ENERGY EC0(MAXGRD)
C
C     -- RHO^13*KX
      DO I=1,2
       FF(I)=KX0(I)
       RK(I)=RHO13(I)*FF(I)
      END DO
C     -- BETA
      SUM_RK=RK(IRA)+RK(IRB)
      PRO_RK=RK(IRA)*RK(IRB)
      BB=QOP*PRO_RK/SUM_RK
C     -- H
      BB2=BB *BB
      BB3=BB2*BB
      BB4=BB3*BB
      HH1=C1*BB+C2
      HH2=BB4 + C3*BB3 + C4*BB2
      HH=HH1/HH2
      TMP0=HH
      RHOAB=RHO(IRA)*RHO(IRB)
      EC0(IIPT)=-RHOAB*TMP0
C
C      1ST DERIVATIVE
C     ------------------
C
C      RHO^13*KX
C
      DO I=1,2
       RHO23(I)=RHO13(I)**2
       DFF(I)=KX(1,I)*X(I)
       DRK_R(I)=CK_R*(FF(I)-DFF(I)*FOUR)/RHO23(I)
       DRK_G(I)=CK_G* DFF(I) *RHO13(I) /GM(I)
      END DO
C     -- BETA
      RKA2=RK(IRA)**2
      RKB2=RK(IRB)**2
      SUM2=SUM_RK**2
      DBDA=QOP*RKB2/SUM2
      DBDB=QOP*RKA2/SUM2
C     -- H
      DHH1=C1
      DHH2=4.0D+00*BB3 + 3.0D+00*C3*BB2 + 2.0D+00*C4*BB
      G1 = DHH1*HH2-HH1*DHH2
      HH22=HH2**2
      DHH= G1/HH22
      DHHDA=DHH*DBDA
      DHHDB=DHH*DBDB
      TMP(IRA)=DHHDA*DRK_R(IRA)
      TMP(IGA)=DHHDA*DRK_G(IRA)
      TMP(IRB)=DHHDB*DRK_R(IRB)
      TMP(IGB)=DHHDB*DRK_G(IRB)
C     **** 1ST DERIVATIVES
      EC(IIPT,IRA)=-RHOAB*TMP(IRA)-RHO(IRB)*TMP0
      EC(IIPT,IRB)=-RHOAB*TMP(IRB)-RHO(IRA)*TMP0
      EC(IIPT,IGA)=-RHOAB*TMP(IGA)
      EC(IIPT,IGB)=-RHOAB*TMP(IGB)
      RHOCHECK=RHOCHECK+RHOTOT*XYZWGT(IPT)
      EOP=EOP+EC0(IIPT)*XYZWGT(IPT)
      IF(NDER.EQ.1) GOTO 100
C
C      2ND DERIVATIVE
C     -------------------
C
C     -- RHO^13*KX
      DO I=1,2
      RHO53(I)=RHO23(I)*RHO(I)
      GM2(I)  =GM(I)**2
      D2FF(I)=KX(2,I)*X(I)*X(I)+KX(1,I)*X(I)
      D2RK_RR(I)=CK_RR*(FF(I)-2.0D0*DFF(I)-8.0D0*D2FF(I))/RHO53(I)
      D2RK_RG(I)=CK_RG*(DFF(I)-4.0D0*D2FF(I))/RHO23(I)/GM(I)
      D2RK_GG(I)=CK_GG*(D2FF(I)-2.0D0*DFF(I))*RHO13(I)/GM2(I)
      END DO
C     -- BETA
      SUM3=SUM2*SUM_RK
      D2BAA=-2.0D0*QOP*RKB2  /SUM3
      D2BBB=-2.0D0*QOP*RKA2  /SUM3
      D2BAB= 2.0D0*QOP*PRO_RK/SUM3
C     -- H
      D2HH2=12.0D0*BB2 + 6.0D0*BB*C3 + 2.0D0*C4
      DG1 =-HH1*D2HH2
      HH23=HH22*HH2
      G2  =DG1*HH2 -2.0D0*G1*DHH2
      D2HH=G2/HH23
      D2HHAA=DHH*D2BAA+D2HH*DBDA*DBDA
      D2HHBB=DHH*D2BBB+D2HH*DBDB*DBDB
      D2HHAB=DHH*D2BAB+D2HH*DBDA*DBDB
      TMP(IRARA)=DHHDA*D2RK_RR(IRA)+D2HHAA*DRK_R(IRA)*DRK_R(IRA)
      TMP(IRAGA)=DHHDA*D2RK_RG(IRA)+D2HHAA*DRK_R(IRA)*DRK_G(IRA)
      TMP(IGAGA)=DHHDA*D2RK_GG(IRA)+D2HHAA*DRK_G(IRA)*DRK_G(IRA)
      TMP(IRBRB)=DHHDB*D2RK_RR(IRB)+D2HHBB*DRK_R(IRB)*DRK_R(IRB)
      TMP(IRBGB)=DHHDB*D2RK_RG(IRB)+D2HHBB*DRK_R(IRB)*DRK_G(IRB)
      TMP(IGBGB)=DHHDB*D2RK_GG(IRB)+D2HHBB*DRK_G(IRB)*DRK_G(IRB)
      TMP(IRARB)=D2HHAB*DRK_R(IRA)*DRK_R(IRB)
      TMP(IGAGB)=D2HHAB*DRK_G(IRA)*DRK_G(IRB)
      TMP(IRAGB)=D2HHAB*DRK_R(IRA)*DRK_G(IRB)
      TMP(IRBGA)=D2HHAB*DRK_G(IRA)*DRK_R(IRB)
C     **** 2ND DERIVATIVES
      EC(IIPT,IRARB)=-RHOAB*TMP(IRARB) -TMP0
     *            -RHO(IRA)*TMP(IRA) -RHO(IRB)*TMP(IRB)
      EC(IIPT,IRARA)=-RHOAB*TMP(IRARA) - RHO(IRB)*TMP(IRA)*2.0D0
      EC(IIPT,IRBRB)=-RHOAB*TMP(IRBRB) - RHO(IRA)*TMP(IRB)*2.0D0
      EC(IIPT,IRAGA)=-RHOAB*TMP(IRAGA) - RHO(IRB)*TMP(IGA)
      EC(IIPT,IRBGB)=-RHOAB*TMP(IRBGB) - RHO(IRA)*TMP(IGB)
      EC(IIPT,IRAGB)=-RHOAB*TMP(IRAGB) - RHO(IRB)*TMP(IGB)
      EC(IIPT,IRBGA)=-RHOAB*TMP(IRBGA) - RHO(IRA)*TMP(IGA)
      EC(IIPT,IGAGA)=-RHOAB*TMP(IGAGA)
      EC(IIPT,IGBGB)=-RHOAB*TMP(IGBGB)
      EC(IIPT,IGAGB)=-RHOAB*TMP(IGAGB)
C
      IF(NDER.EQ.2) GOTO 100
C
C      3RD DERIVATIVE
C     ----------------
C
      DO I=1,2
C      RHO^13*KX
      RHO83(I)=RHO53(I)*RHO(I)
      GM3(I)  =GM2(I)*GM(I)
      D3FF(I) =KX(3,I)*X(I)**3 + THREE*KX(2,I)*X(I)**2 + KX(1,I)*X(I)
C
      D3RK_RRR(I)=CK_RRR*(FIVE*FF(I)-SIX*DFF(I)-SIX*EIGHT*D2FF(I)
     * -FOUR*EIGHT*D3FF(I))/RHO83(I)
      D3RK_RRG(I)=CK_RRG*(DFF(I) -TWO*D2FF(I) -EIGHT*D3FF(I))/GM(I)
     * /RHO53(I)
      D3RK_RGG(I)=CK_RGG*(TWO*DFF(I)-XNINE*D2FF(I)+FOUR*D3FF(I))/GM2(I)
     * /RHO23(I)
      D3RK_GGG(I)=CK_GGG*(EIGHT*DFF(I)-SIX*D2FF(I)+D3FF(I))/GM3(I)
     * *RHO13(I)
      END DO
C      BETA
      SUM4=SUM3*SUM_RK
      D3BAAA=SIX*QOP*RKB2/SUM4
      D3BBBB=SIX*QOP*RKA2/SUM4
      D3BAAB=TWO*QOP*(RKB2-TWO*PRO_RK)/SUM4
      D3BABB=TWO*QOP*(RKA2-TWO*PRO_RK)/SUM4
C      H
      D3HH2=SIX*FOUR*BB + SIX*C3
      D2G1= -DHH1*D2HH2 - HH1*D3HH2
      DG2 = D2G1*HH2 - DG1*DHH2 - TWO*G1*D2HH2
      HH24=HH22**2
      D3HH= (DG2*HH2 -THREE*G2*DHH2)/HH24
      D3HHAAA=D3HH*DBDA**3 + THREE*D2HH*D2BAA*DBDA + DHH*D3BAAA
      D3HHBBB=D3HH*DBDB**3 + THREE*D2HH*D2BBB*DBDB + DHH*D3BBBB
      D3HHAAB=D3HH*(DBDA**2)*DBDB + DHH*D3BAAB
     *      +D2HH*(TWO*D2BAB*DBDA + D2BAA*DBDB)
      D3HHABB=D3HH*(DBDB**2)*DBDA + DHH*D3BABB
     *      +D2HH*(TWO*D2BAB*DBDB + D2BBB*DBDA)
C
      TMP(IRARARA)=DHHDA*D3RK_RRR(IRA)+D3HHAAA*DRK_R(IRA)*DRK_R(IRA)*
     * DRK_R(IRA)+THREE*D2HHAA*D2RK_RR(IRA)*DRK_R(IRA)
      TMP(IGAGAGA)=DHHDA*D3RK_GGG(IRA)+D3HHAAA*DRK_G(IRA)*DRK_G(IRA)*
     * DRK_G(IRA)+THREE*D2HHAA*D2RK_GG(IRA)*DRK_G(IRA)
      TMP(IRARAGA)=DHHDA*D3RK_RRG(IRA)+D3HHAAA*DRK_R(IRA)*DRK_R(IRA)*
     * DRK_G(IRA)+D2HHAA*(D2RK_RR(IRA)*DRK_G(IRA)+TWO*D2RK_RG(IRA)*
     * DRK_R(IRA))
      TMP(IRAGAGA)=DHHDA*D3RK_RGG(IRA)+D3HHAAA*DRK_R(IRA)*DRK_G(IRA)*
     * DRK_G(IRA)+D2HHAA*(D2RK_GG(IRA)*DRK_R(IRA)+TWO*D2RK_RG(IRA)*
     * DRK_G(IRA))
C
      TMP(IRBRBRB)=DHHDB*D3RK_RRR(IRB)+D3HHBBB*DRK_R(IRB)*DRK_R(IRB)*
     * DRK_R(IRB)+THREE*D2HHBB*D2RK_RR(IRB)*DRK_R(IRB)
      TMP(IGBGBGB)=DHHDB*D3RK_GGG(IRB)+D3HHBBB*DRK_G(IRB)*DRK_G(IRB)*
     * DRK_G(IRB)+THREE*D2HHBB*D2RK_GG(IRB)*DRK_G(IRB)
      TMP(IRBRBGB)=DHHDB*D3RK_RRG(IRB)+D3HHBBB*DRK_R(IRB)*DRK_R(IRB)*
     * DRK_G(IRB)+D2HHBB*(D2RK_RR(IRB)*DRK_G(IRB)+TWO*D2RK_RG(IRB)*
     * DRK_R(IRB))
      TMP(IRBGBGB)=DHHDB*D3RK_RGG(IRB)+D3HHBBB*DRK_R(IRB)*DRK_G(IRB)*
     * DRK_G(IRB)+D2HHBB*(D2RK_GG(IRB)*DRK_R(IRB)+TWO*D2RK_RG(IRB)*
     * DRK_G(IRB))
C
      TMP(IRARARB)=D3HHAAB*DRK_R(IRA)*DRK_R(IRA)*DRK_R(IRB)
     *          +D2HHAB*D2RK_RR(IRA)*DRK_R(IRB)
      TMP(IRARAGB)=D3HHAAB*DRK_R(IRA)*DRK_R(IRA)*DRK_G(IRB)
     *          +D2HHAB*D2RK_RR(IRA)*DRK_G(IRB)
      TMP(IRARBGA)=D3HHAAB*DRK_R(IRA)*DRK_G(IRA)*DRK_R(IRB)
     *          +D2HHAB*D2RK_RG(IRA)*DRK_R(IRB)
      TMP(IRAGAGB)=D3HHAAB*DRK_R(IRA)*DRK_G(IRA)*DRK_G(IRB)
     *          +D2HHAB*D2RK_RG(IRA)*DRK_G(IRB)
      TMP(IRBGAGA)=D3HHAAB*DRK_G(IRA)*DRK_G(IRA)*DRK_R(IRB)
     *          +D2HHAB*D2RK_GG(IRA)*DRK_R(IRB)
      TMP(IGAGAGB)=D3HHAAB*DRK_G(IRA)*DRK_G(IRA)*DRK_G(IRB)
     *          +D2HHAB*D2RK_GG(IRA)*DRK_G(IRB)
C
      TMP(IRARBRB)=D3HHABB*DRK_R(IRB)*DRK_R(IRB)*DRK_R(IRA)
     *          +D2HHAB*D2RK_RR(IRB)*DRK_R(IRA)
      TMP(IRBRBGA)=D3HHABB*DRK_R(IRB)*DRK_R(IRB)*DRK_G(IRA)
     *          +D2HHAB*D2RK_RR(IRB)*DRK_G(IRA)
      TMP(IRARBGB)=D3HHABB*DRK_R(IRB)*DRK_G(IRB)*DRK_R(IRA)
     *          +D2HHAB*D2RK_RG(IRB)*DRK_R(IRA)
      TMP(IRBGAGB)=D3HHABB*DRK_R(IRB)*DRK_G(IRB)*DRK_G(IRA)
     *          +D2HHAB*D2RK_RG(IRB)*DRK_G(IRA)
      TMP(IRAGBGB)=D3HHABB*DRK_G(IRB)*DRK_G(IRB)*DRK_R(IRA)
     *          +D2HHAB*D2RK_GG(IRB)*DRK_R(IRA)
      TMP(IGAGBGB)=D3HHABB*DRK_G(IRB)*DRK_G(IRB)*DRK_G(IRA)
     *          +D2HHAB*D2RK_GG(IRB)*DRK_G(IRA)
C
      EC(IIPT,IRARARB)=EC(IIPT,IRARARB)-RHOAB*TMP(IRARARB)-TMP(IRA)*TWO
     *            -RHO(IRA)*TMP(IRARA) - RHO(IRB)*TMP(IRARB)*TWO
      EC(IIPT,IRARBRB)=EC(IIPT,IRARBRB)-RHOAB*TMP(IRARBRB)-TMP(IRB)*TWO
     *            -RHO(IRB)*TMP(IRBRB) - RHO(IRA)*TMP(IRARB)*TWO
      EC(IIPT,IRARBGA)=EC(IIPT,IRARBGA)-RHOAB*TMP(IRARBGA) - TMP(IGA)
     *            -RHO(IRA)*TMP(IRAGA) - RHO(IRB)*TMP(IRBGA)
      EC(IIPT,IRARBGB)=EC(IIPT,IRARBGB)-RHOAB*TMP(IRARBGB) - TMP(IGB)
     *            -RHO(IRA)*TMP(IRAGB) - RHO(IRB)*TMP(IRBGB)
C
      EC(IIPT,IRARARA)=EC(IIPT,IRARARA)-RHOAB*TMP(IRARARA) -
     *  RHO(IRB)*TMP(IRARA)*THREE
      EC(IIPT,IRBRBRB)=EC(IIPT,IRBRBRB)-RHOAB*TMP(IRBRBRB) -
     *  RHO(IRA)*TMP(IRBRB)*THREE
C
      EC(IIPT,IRARAGA)=EC(IIPT,IRARAGA)-RHOAB*TMP(IRARAGA) -
     * RHO(IRB)*TMP(IRAGA)*TWO
      EC(IIPT,IRBRBGB)=EC(IIPT,IRBRBGB)-RHOAB*TMP(IRBRBGB) -
     * RHO(IRA)*TMP(IRBGB)*TWO
      EC(IIPT,IRARAGB)=EC(IIPT,IRARAGB)-RHOAB*TMP(IRARAGB) -
     * RHO(IRB)*TMP(IRAGB)*TWO
      EC(IIPT,IRBRBGA)=EC(IIPT,IRBRBGA)-RHOAB*TMP(IRBRBGA) -
     * RHO(IRA)*TMP(IRBGA)*TWO
C
      EC(IIPT,IRAGAGA)=EC(IIPT,IRAGAGA)-RHOAB*TMP(IRAGAGA)
     * -RHO(IRB)*TMP(IGAGA)
      EC(IIPT,IRBGBGB)=EC(IIPT,IRBGBGB)-RHOAB*TMP(IRBGBGB)
     * -RHO(IRA)*TMP(IGBGB)
      EC(IIPT,IRAGAGB)=EC(IIPT,IRAGAGB)-RHOAB*TMP(IRAGAGB)
     * -RHO(IRB)*TMP(IGAGB)
      EC(IIPT,IRBGAGA)=EC(IIPT,IRBGAGA)-RHOAB*TMP(IRBGAGA)
     * -RHO(IRA)*TMP(IGAGA)
      EC(IIPT,IRBGAGB)=EC(IIPT,IRBGAGB)-RHOAB*TMP(IRBGAGB)
     * -RHO(IRA)*TMP(IGAGB)
      EC(IIPT,IRAGBGB)=EC(IIPT,IRAGBGB)-RHOAB*TMP(IRAGBGB)
     * -RHO(IRB)*TMP(IGBGB)
C
      EC(IIPT,IGAGAGA)=EC(IIPT,IGAGAGA)-RHOAB*TMP(IGAGAGA)
      EC(IIPT,IGBGBGB)=EC(IIPT,IGBGBGB)-RHOAB*TMP(IGBGBGB)
      EC(IIPT,IGAGAGB)=EC(IIPT,IGAGAGB)-RHOAB*TMP(IGAGAGB)
      EC(IIPT,IGAGBGB)=EC(IIPT,IGAGBGB)-RHOAB*TMP(IGAGBGB)
C     -- DENSITY CUT-OFF
      ENDIF
  100 CONTINUE
      END DO
      ECOR=ECOR+EOP
C
      IF(DBUG) THEN
      WRITE(6,*)'INTEGRATED RHO IN DEROPD=',RHOCHECK
      WRITE(6,*)'OP CORRELATION ENERGY=',EOP
      ENDIF
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFVWN
      SUBROUTINE NDFVWN(ECOR,EC0,EC,DRHOI_,XYZWGT,SCALE,RHOMIN,
     * ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL DBUG
      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 DRHOI_(NPTGRD,4,IRU)
      DIMENSION EC0(*)
      DIMENSION EC(ILENG,*)
      DIMENSION XYZWGT(*)
      DIMENSION A(3),B(3),C(3),X0(3),Q(3)
      DIMENSION X0I(3),C1(3),C2(3),C3(3),XI(3),XXI(3),
     *          DI0(3),DI1(3),DI12(3),DI4(3),DENOM(3)
      DIMENSION E(3),EDX(3),ED(3),EDDX(3),EDD(3),E3DX(3),E3D(3)
      DIMENSION T1D(3),T2D(3),T3D(3)
      DATA ZERO /0.0D+00/
      DATA PT5  /0.5D+00/
      DATA ONE  /1.0D+00/
      DATA TWO  /2.0D+00/
      DATA THREE/3.0D+00/
      DATA FOUR /4.0D+00/
      DATA FIVE /5.0D+00/
      DATA SIX  /6.0D+00/
      DATA SEVEN/7.0D+00/
      DATA ANINE/9.0D+00/
      DBUG=.FALSE.
C
      SIXTH=ONE/SIX
C      VARIVABLES AND CONSTANTS FOR VWN5 FUNCITONAL
C     ---------------------------------------------
      FACR = 0.78762331789974325053D+00
      FACH = 1.70992093416136561756D+00
      GAMMA= 0.5198420997897464D+00
C       VWN'S CONSTANT PARAMETERS
C      --------------------------
      VWNAP=0.0621814D+00/TWO
      VWNBP=3.7274400D+00
      VWNCP=12.9352000D+00
      VWNX0P= -0.1049800D+00
      VWNAF=VWNAP/TWO
      VWNBF=7.0604200D+00
      VWNCF=18.0578000D+00
      VWNX0F=-0.3250000D+00
      VWNAA=-0.03377372D+00/TWO
      VWNBA= 1.1310700D+00
      VWNCA= 13.0045000D+00
      VWNX0A= -0.0047584D+00
C
      EVWN=ZERO
      RHOCHECK=ZERO
      DO IPT=IST,IEND
      IIPT=IPT-IST+1
      RHOA=DRHOI_(IPT,1,1)
      IF(IRU.EQ.1) THEN
       RHOB=RHOA
      ELSE
       RHOB=DRHOI_(IPT,1,IRU)
      ENDIF
      RHOT = RHOA + RHOB
      IF (RHOT/TWO.GT.RHOMIN) THEN
C           RHOB  = RO(IGRD,RB)
C           RHOT  = RHOA + RHOB
           RHO16 = RHOT**SIXTH
           A(1) = VWNAP
           A(2) = VWNAF
           A(3) = VWNAA
           B(1) = VWNBP
           B(2) = VWNBF
           B(3) = VWNBA
           C(1) = VWNCP
           C(2) = VWNCF
           C(3) = VWNCA
           X0(1)= VWNX0P
           X0(2)= VWNX0F
           X0(3)= VWNX0A
           Q(1) = 6.15199081975908D+00
           Q(2) = 4.730926909560114D+00
           Q(3) = 7.123108759029235D+00
           ZETA  = (RHOA - RHOB)/RHOT
           ZETA2 = ZETA *ZETA
           ZETA3 = ZETA2*ZETA
           ZETA4 = ZETA2*ZETA2
           ONEPZ= ONE + ZETA
           ONEMZ= ONE - ZETA
           OPZPP= ONEPZ**(-TWO/THREE)
           OPZP = OPZPP*ONEPZ
           OPZ  = OPZP*ONEPZ
           OMZPP= ONEMZ**(-TWO/THREE)
           IF(RHOB.EQ.ZERO) OMZPP=ZERO
           OMZP = OMZPP*ONEMZ
           OMZ  = OMZP*ONEMZ
           FT2  = OPZ + OMZ
           F    = (FT2 - TWO)/GAMMA
           Y    = ZETA4
           V    = (ONE - ZETA4)/FACH
           X    = FACR/RHO16
           X2   = X*X
           DO I = 1, 3
              X0I(I)= X0(I)*X0(I)+B(I)*X0(I)+C(I)
              C1(I)= TWO*B(I)/Q(I)
              C2(I)= B(I)*X0(I)/X0I(I)
              C3(I)= C1(I)+FOUR*X0(I)/Q(I)
              XI(I)= X2+B(I)*X+C(I)
              AT1= LOG(X2/XI(I))
              DI1(I)= TWO*X+B(I)
              AT2= DATAN(Q(I)/DI1(I))
              AT3= LOG((X-X0(I))*(X-X0(I))/XI(I))
              E(I)= A(I)*(AT1+C1(I)*AT2-C2(I)*(AT3+C3(I)*AT2))
           END DO
           D     = E(2) - E(1)
           TDEC0 = (E(3)*V + Y*D)
           ADEC  = F*TDEC0
           ECF00 = E(1)+ADEC
C
C           VWN ENERGY
C           -----------
           EC0(IIPT) = EC0(IIPT)+SCALE*ECF00*RHOT
C
C            1ST DERIVATIVE
C            --------------
           P2   = RHOT*RHOT
           ZA   = TWO*RHOB/P2
           ZB   = -TWO*RHOA/P2
           GAM43= FOUR/(THREE*GAMMA)
           FT1  = OPZP - OMZP
           FT3  = GAM43*FT1
           FA   = ZA*FT3
           FB   = ZB*FT3
C
           DXDP = -SIXTH*X/RHOT
           DO I = 1, 3
              XXI(I)=X*XI(I)
              DI0(I)=TWO*C(I)+B(I)*X
              T1D(I)=DI0(I)/XXI(I)
              DI12(I)=DI1(I)*DI1(I)
              DENOM(I)=Q(I)*Q(I)+DI12(I)
              T2D(I)=-TWO*Q(I)/DENOM(I)
              DI4(I)=(X-X0(I))*XI(I)
              T3D(I)=(DI0(I)+DI1(I)*X0(I))/DI4(I)
              EDX(I)=A(I)*(T1D(I)+C1(I)*T2D(I)-C2(I)*
     *       (T3D(I)+C3(I)*T2D(I)))
              ED(I)=DXDP*EDX(I)
           END DO
C
           DD = ED(2) - ED(1)
           YA = FOUR*ZETA3*ZA
           YB = FOUR*ZETA3*ZB
           VA = -YA/FACH
           VB = -YB/FACH
C
           TDEC6 = ED(3)*V + Y*DD
           TDEC1 = (E(3)*VA + YA*D+TDEC6)
           TDEC2 = (E(3)*VB + YB*D+TDEC6)
           DECA  = FA*TDEC0 + F*TDEC1
           DECB  = FB*TDEC0 + F*TDEC2
C
C          ASSEMBLE 1ST DERIVATIVES
C            ----------------------------
           ECF10 = DECA + ED(1)
           ECF01 = DECB + ED(1)
C
C            1ST DERIVATIVES OF THE VWN FUNCTIONAL
C            ----------------------------------------
           EC(IIPT,IRA) = EC(IIPT,IRA) + SCALE*(ECF00+RHOT*ECF10)
           EC(IIPT,IRB) = EC(IIPT,IRB) + SCALE*(ECF00+RHOT*ECF01)
C     -- DEBUG DENSITY
           RHOCHECK=RHOCHECK+RHOT*XYZWGT(IPT)
           EVWN=EVWN+SCALE*ECF00*RHOT*XYZWGT(IPT)
           IF(NDER.EQ.1) GOTO 100
C
C            2ND DERIVATIVE
C            --------------
           P3   = P2*RHOT
           GAM49= FOUR/(ANINE*GAMMA)
           ZA2  = ZA*ZA
           ZB2  = ZB*ZB
           ZAB  = TWO*(RHOA-RHOB)/P3
           ZAA  = -FOUR*RHOB/P3
           ZBB  = FOUR*RHOA/P3
C
           FT0  = OPZPP + OMZPP
           FT4  = GAM49*FT0
           FAA  = ZAA*FT3 + FT4*ZA2
           FBB  = ZBB*FT3 + FT4*ZB2
           FAB  = ZAB*FT3 + FT4*ZA*ZB
C
           D2XDP2 = -SEVEN*SIXTH*DXDP/RHOT
           DXDP2=DXDP*DXDP
           DO I = 1, 3
              T1DD=(-T1D(I)*(C(I)+TWO*B(I)*X+THREE*X2)+B(I))/XXI(I)
              T2DD=-T2D(I)*FOUR*DI1(I)/DENOM(I)
              T3DD=(B(I)+TWO*X0(I))/DI4(I)-T3D(I)/(X-X0(I))-DI1(I)
     *        *T3D(I)/XI(I)
              EDDX(I)=A(I)*(T1DD+C1(I)*T2DD-C2(I)*(T3DD+C3(I)*T2DD))
              EDD(I)=DXDP2*EDDX(I)+D2XDP2*EDX(I)
           END DO
C
C            DERIVATIVES OF DELTAEC
C            ----------------------
           YAA  = FOUR*ZETA2*(THREE*ZA2 + ZETA*ZAA)
           YBB  = FOUR*ZETA2*(THREE*ZB2 + ZETA*ZBB)
           YAB  = FOUR*ZETA2*(THREE*ZA*ZB + ZETA*ZAB)
           DDD  = EDD(2)- EDD(1)
           VAB  = -YAB/FACH
           VAA  = -YAA/FACH
           VBB  = -YBB/FACH
           TDEC3 = EDD(3)*V+Y*DDD
           TDEC4 = TWO*(ED(3)*VA + YA*DD)
           TDEC5 = TWO*(ED(3)*VB + YB*DD)
           DECAA = FAA*TDEC0+TWO*FA*TDEC1+F*(TDEC4+TDEC3+E(3)*VAA+YAA*D)
           DECBB = FBB*TDEC0+TWO*FB*TDEC2+F*(TDEC5+TDEC3+E(3)*VBB+YBB*D)
           DECAB = FAB*TDEC0+FA*TDEC2+FB*TDEC1+F*(TDEC3+PT5*
     *     (TDEC4+TDEC5)+E(3)*VAB+YAB*D)
C
C            FINALLY ASSEMBLE DERIVATIVES
C            ----------------------------
           ECF20 = DECAA+ EDD(1)
           ECF02 = DECBB+ EDD(1)
           ECF11 = DECAB+ EDD(1)
C
C            SECOND DERIVATIVES OF THE VWN FUNCTIONAL
C            ----------------------------------------
           EC(IIPT,IRARA) =EC(IIPT,IRARA)+ SCALE*(TWO*ECF10+RHOT*ECF20)
           EC(IIPT,IRARB) =EC(IIPT,IRARB)+ SCALE*
     *     (ECF10 + ECF01 + RHOT*ECF11)
           EC(IIPT,IRBRB) =EC(IIPT,IRBRB)+
     *     SCALE*(TWO*ECF01 + RHOT*ECF02)
           IF(NDER.EQ.2) GOTO 100
C
C            THIRD DERIVATIVE
C            -----------------
           ZA3  = ZA2*ZA
           ZB3  = ZB2*ZB
           P4   = P3*RHOT
           ZAAA = THREE*FOUR*RHOB/P4
           ZAAB = FOUR*(TWO*RHOB-RHOA)/P4
           ZABB = -FOUR*(TWO*RHOA-RHOB)/P4
           ZBBB = -THREE*FOUR*RHOA/P4
           FT5 = GAM49 *(-TWO/THREE)*((1+ZETA)**(-FIVE/THREE)
     *     -(1-ZETA)**(-FIVE/THREE))
           FAAA = FT5*ZA3+THREE*FT4*ZA*ZAA+FT3*ZAAA
           FAAB = FT5*ZA2*ZB+FT4*(TWO*ZA*ZAB+ZAA*ZB)+FT3*ZAAB
           FABB = FT5*ZA*ZB2+FT4*(TWO*ZB*ZAB+ZA*ZBB)+FT3*ZABB
           FBBB = FT5*ZB3+THREE*FT4*ZB*ZBB+FT3*ZBBB
           D3XDP3 = -13.0D+00*SIXTH*D2XDP2/RHOT
           X3     = X2*X
           DXDP3=DXDP2*DXDP
           DO I=1,3
              DI3= B(I)+TWO*X0(I)
              DENOM2=DENOM(I)*DENOM(I)
              XI2=XI(I)*XI(I)
              XI3=XI2*XI(I)
              T1DDD= (-TWO*B(I)/X2+TWO*DI0(I)/X3)/XI(I)-(TWO*DI1(I)
     *        *(B(I)/X-DI0(I)/X2))/XI2+(DI0(I)
     *        *(TWO*DI12(I)/XI3-TWO/XI2))/X
              T2DDD=FOUR*FOUR/DENOM2*(-FOUR*Q(I)*DI12(I)/DENOM(I)+Q(I))
              T3DDD=-TWO/(X-X0(I))/XI(I)*(DI3/(X-X0(I))+DI1(I)*
     *        DI3/XI(I)-(DI0(I)+DI1(I)*X0(I))/(X-X0(I))/(X-X0(I))
     *       -DI1(I)*(DI0(I)+DI1(I)*X0(I))/XI(I)/(X-X0(I))-DI12(I)*
     *        (DI0(I)+DI1(I)*X0(I))/XI2+(DI0(I)+DI1(I)*X0(I))/XI(I))
              E3DX(I)=A(I)*(T1DDD+C1(I)*T2DDD-C2(I)*(T3DDD+C3(I)*T2DDD))
              E3D(I)=EDX(I)*D3XDP3+THREE*EDDX(I)*DXDP*D2XDP2+E3DX(I)
     *        *DXDP3
           END DO
           D3D   = E3D(2)-E3D(1)
           YAAA  = FOUR*(SIX*ZETA*ZA3+ANINE*ZETA2*ZA*ZAA+ZETA3*ZAAA)
           YAAB  = FOUR*(SIX*ZETA*ZA2*ZB+THREE*ZETA2*
     *     (TWO*ZA*ZAB+ZB*ZAA)+ZETA3*ZAAB)
           YABB  = FOUR*(SIX*ZETA*ZA*ZB2+THREE*ZETA2*
     *     (ZA*ZBB+TWO*ZB*ZAB)+ZETA3*ZABB)
           YBBB  = FOUR*(SIX*ZETA*ZB3+ANINE*ZETA2*ZB*ZBB+ZETA3*ZBBB)
           VAAA  = -YAAA/FACH
           VAAB  = -YAAB/FACH
           VABB  = -YABB/FACH
           VBBB  = -YBBB/FACH
           TDEC7 = V*E3D(3)+Y*D3D
           TDEC8 = ED(3)*VAA+DD*YAA
           TDEC9 = ED(3)*VBB+DD*YBB
           TDEC10= ED(3)*VAB+DD*YAB
           DECAAA = FAAA*TDEC0+THREE*FAA*TDEC1+THREE*FA*
     *     (TDEC4+TDEC3+E(3)*VAA+YAA*D)+F*(E(3)*VAAA+D*YAAA+
     *     THREE*TDEC8+THREE*(EDD(3)*VA+DDD*YA)+TDEC7)
           DECAAB = FAAB*TDEC0+TWO*FAB*TDEC1+FAA*TDEC2+FB*
     *     (TDEC4+TDEC3+E(3)*VAA+D*YAA)+TWO*FA*(TDEC3+PT5*
     *     (TDEC4+TDEC5)+E(3)*VAB+D*YAB)+F*(EDD(3)*VB+DDD*YB
     *     +TWO*TDEC10+TDEC7+TDEC8+TWO*(EDD(3)*VA+DDD*YA)
     *     +E(3)*VAAB+D*YAAB)
           DECABB = FABB*TDEC0+TWO*FAB*TDEC2+FBB*TDEC1+FA*(TDEC3+
     *     TDEC5+E(3)*VBB+D*YBB)+TWO*FB*(TDEC3+PT5*(TDEC4+TDEC5)
     *     +E(3)*VAB+D*YAB)+F*(EDD(3)*VA+DDD*YA+TWO*TDEC10+TDEC7
     *     +TDEC9+TWO*(EDD(3)*VB+DDD*YB)+E(3)*VABB+D*YABB)
           DECBBB = FBBB*TDEC0+THREE*FBB*TDEC2+THREE*FB*(TDEC5+
     *     TDEC3+E(3)*VBB+YBB*D)+F*(E(3)*VBBB+D*YBBB+THREE*TDEC9
     *     +THREE*(EDD(3)*VB+DDD*YB)+TDEC7)
C
C            FINALLY ASSEMBLE DERIVATIVES
C            ----------------------------
           ECF30 = DECAAA+E3D(1)
           ECF21 = DECAAB+E3D(1)
           ECF12 = DECABB+E3D(1)
           ECF03 = DECBBB+E3D(1)
C
C            THIRD DERIVATIVES OF THE VWN FUNCTIONAL
C            ----------------------------------------
           EC(IIPT,IRARARA)=EC(IIPT,IRARARA)+SCALE*
     *     (THREE*ECF20+ECF30*RHOT)
           EC(IIPT,IRARARB)=EC(IIPT,IRARARB)+SCALE*
     *     (TWO*ECF11+ECF20+ECF21*RHOT)
           EC(IIPT,IRARBRB)=EC(IIPT,IRARBRB)+SCALE*
     *     (TWO*ECF11+ECF02+ECF12*RHOT)
           EC(IIPT,IRBRBRB)=EC(IIPT,IRBRBRB)+SCALE*
     *     (THREE*ECF02+ECF03*RHOT)
      ENDIF
  100 CONTINUE
      END DO
      ECOR=ECOR+EVWN
C
      IF(DBUG) THEN
      WRITE(6,*)'INTEGRATED RHO IN DERVWN=',RHOCHECK
      WRITE(6,*)'VWN CORRELATION ENERGY=',EVWN
      ENDIF
C
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFLYP
      SUBROUTINE NDFLYP(ECOR,EC0,EC,DRHOI_,XYZWGT,SCALE,RHOMIN,
     * ILENG,NPTGRD,NDER,IRU,IST,IEND)
C       _________________________________________________
C
C
C       LEE-YANG-PARR (LYP) CORRELATION FUNCTIONAL
C
C            C. LEE, W. YANG, AND R. G. PARR
C            PHYS. REV. B 37, 785-789 (1988)
C
C       __________________________________________________
C               THE NECESSARY FORMULAE ARE SPELLED OUT IN
C                        B.G.JOHNSON, P.M.W.GILL, J.A.POPLE
C                            J.CHEM.PHYS. 98, 5612-26(1993)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION MBABOAAA,MBABOABB,MBABOAAB,MBABOBBB,MBABOBAA,
     * MBABOBAB,MBABOAAADA,MBABOABBDA,MBABOAABDA,MBABOBBBDA,MBABOBAADA,
     * MBABOBABDA,MBABOAAADB,MBABOABBDB,MBABOAABDB,MBABOBBBDB,MBABOBAADB
     * ,MBABOBABDB,MBTERM2A,MBTERM2B,MBTERM2ADA, MBTERM2ADB,MBTERM2BDA,
     *  MBTERM2BDB
      LOGICAL DBUG
      PARAMETER (ALYPA=0.04918D+00,ALYPB=0.13200D+00,ALYPC=0.25330D+00,
     * ALYPD=0.34900D+00,C1=1D+00,C2=2D+00,C3=3D+00,C4=4D+00,C6=6D+00,
     * C7=7D+00,C8=8.0D+00,C9=9.0D+00,C11=11D+00,C47=47D+00,C40=40D+00,
     * C28=28D+00, C5 = 5D+00)
      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 EC(ILENG,*)
      DIMENSION EC0(*)
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION XYZWGT(*)
      DATA ZERO/0.0D+00/
      DATA TWO /2.0D+00/
      DBUG=.FALSE.
C
C     CONSTANTS
C     ---------------
      PI=ATAN(1.0D+00)*4.0D+00
      FACE  = C2**(C11/C3) * 0.3D+00 * (C3*PI**2)**(C2/C3)
      ALYPAB= ALYPA * ALYPB
      ALYPD2 = ALYPD*ALYPD
C     DENSITY SETTING
C     ---------------
      RHOCHECK=ZERO
      E_LYP=ZERO
      DO IPT=IST,IEND
      IIPT=IPT-IST+1
      RHOA  =DRHOI_(IPT,1,1)
      RHOB=DRHOI_(IPT,1,IRU)
      RHOT  =RHOA+RHOB
      IF(RHOT/TWO.GT.RHOMIN) THEN
      DXA=DRHOI_(IPT,2,1)
      DYA=DRHOI_(IPT,3,1)
      DZA=DRHOI_(IPT,4,1)
      DXB=DRHOI_(IPT,2,IRU)
      DYB=DRHOI_(IPT,3,IRU)
      DZB=DRHOI_(IPT,4,IRU)
      GRDAA=DXA*DXA + DYA*DYA + DZA*DZA
      GRDBB=DXB*DXB + DYB*DYB + DZB*DZB
      GRDAB=DXA*DXB + DYA*DYB + DZA*DZB
      RHO13 =RHOT**(1.0D+00/C3)
      RHO13A=RHOA**(1.0D+00/C3)
      RHO13B=RHOB**(1.0D+00/C3)
      RHO23 = RHO13*RHO13
      RHO2  =RHOT * RHOT
      RHO2A =RHOA * RHOA
      RHO2B =RHOB * RHOB
      RHO3 = RHO2*RHOT
      RHO43 =RHO13 * RHOT
      RHO43A=RHO13A* RHOA
      RHO43B=RHO13B* RHOB
      RHO53 = RHO43*RHO13
      RHO73 = RHO53*RHO23
      RHO83=RHO43* RHO43
      RHO83A=RHO43A* RHO43A
      RHO83B=RHO43B* RHO43B
      RHO113=RHO43 * RHO43 * RHOT
      RHO133=RHO113*RHO23
      RHO103=RHO73*RHOT
C
C      DELTA AND OMEGA
C      ---------------
      XC    =ALYPC/RHO13
      XD    =ALYPD/RHO13
      XDNM  =1.0D+00+XD
      XDNM2 = XDNM*XDNM
      RX    =XD/XDNM
      DELTA =XC+RX
      OMEGA =DEXP(-XC)/XDNM/RHO113
      ABOMEGA=ALYPAB*OMEGA
C
C      LYP ENERGY TERM
C      ---------------
      TERM1A=(1.0D+00-C3*DELTA-(DELTA-C11)*RHOA/RHOT)/C9
      TERM1B=(1.0D+00-C3*DELTA-(DELTA-C11)*RHOB/RHOT)/C9
      DFDAA=-ABOMEGA*(RHOA*RHOB*TERM1A-RHO2B)
      DFDBB=-ABOMEGA*(RHOA*RHOB*TERM1B-RHO2A)
      DFDAB=-ABOMEGA*(RHOA*RHOB*(C47-C7*DELTA)/C9-C4/C3*RHO2)
      TERM2 = -C4*ALYPA/XDNM*RHOA*RHOB/RHOT
      TERM3 = RHOA*RHOB*(RHO83A+RHO83B)
      ELYP    =  TERM2-FACE*ABOMEGA*TERM3+DFDAA*GRDAA+DFDBB*GRDBB+
     * DFDAB*GRDAB
      EC0(IIPT) =EC0(IIPT)+ ELYP*SCALE
C
C      1ST DERIVATIVE
C      --------------
      DOMEGA1 =-(C11*RHO13-ALYPC-ALYPD/XDNM)/RHO43/C3
      DOMEGA = DOMEGA1*OMEGA
      DDELTA=(ALYPD*ALYPD /RHOT/RHO13/RHO13 /XDNM/XDNM - DELTA/RHOT) /C3
      TERM4A=(C3+RHOA/RHOT)*DDELTA
      TERM4B=(C3+RHOB/RHOT)*DDELTA
      TERM5A =(DELTA-C11) * RHOA/RHO2
      TERM5B =(DELTA-C11) * RHOB/RHO2
      TERM6A=C11/C3*RHO83A+RHO83B
      TERM6B=C11/C3*RHO83B+RHO83A
C
C      PREPARE FOR THE 2ND DERIVATIVES
C      -------------------------------
      MBABOAAA=RHOB*TERM1A-RHOA*RHOB/C9*(TERM4A + TERM5B)
      MBABOABB=RHOB*TERM1B-RHOA*RHOB/C9*(TERM4B - TERM5B)-C2*RHOA
      MBABOAAB=RHOB*(C47-C7*DELTA)/C9 - C7/C9*RHOA*RHOB*DDELTA -
     * C8/C3*RHOT
      MBABOBBB=RHOA*TERM1B-RHOA*RHOB/C9*(TERM4B + TERM5A)
      MBABOBAA=RHOA*TERM1A-RHOA*RHOB/C9*(TERM4A - TERM5A)-C2*RHOB
      MBABOBAB=RHOA*(C47-C7*DELTA)/C9 - C7/C9*RHOA*RHOB*DDELTA -
     * C8/C3*RHOT
      MBTERM2A = (ALYPD/RHO43/XDNM/C3 + 1.0D+00/RHOA-1.0D+00/RHOT)
      MBTERM2B = (ALYPD/RHO43/XDNM/C3 + 1.0D+00/RHOB-1.0D+00/RHOT)
C      --------------------------------
      D2AFDAA=DOMEGA1*DFDAA - ABOMEGA*MBABOAAA
      D2AFDBB=DOMEGA1*DFDBB - ABOMEGA*MBABOABB
      D2AFDAB=DOMEGA1*DFDAB - ABOMEGA*MBABOAAB
      D2BFDBB=DOMEGA1*DFDBB - ABOMEGA*MBABOBBB
      D2BFDAA=DOMEGA1*DFDAA - ABOMEGA*MBABOBAA
      D2BFDAB=DOMEGA1*DFDAB - ABOMEGA*MBABOBAB
C      ---------------
      DEDRA  =TERM2 * MBTERM2A
     *       -FACE*ALYPAB*(DOMEGA*TERM3+OMEGA*RHOB*TERM6A)
     *       +D2AFDAA*GRDAA + D2AFDBB*GRDBB + D2AFDAB*GRDAB
      DEDRB  =TERM2 * MBTERM2B
     *       -FACE*ALYPAB*(DOMEGA*TERM3+OMEGA*RHOA*TERM6B)
     *       +D2BFDAA*GRDAA + D2BFDBB*GRDBB + D2BFDAB*GRDAB
C
C      SCALE PARAMETER FOR HYBRID FUNCTIONAL
C      -------------------------------------
      EC(IIPT,IRA) = EC(IIPT,IRA) + SCALE * DEDRA
      EC(IIPT,IRB) = EC(IIPT,IRB) + SCALE * DEDRB
      EC(IIPT,IGA) = EC(IIPT,IGA) + SCALE * DFDAA
      EC(IIPT,IGB) = EC(IIPT,IGB) + SCALE * DFDBB
      EC(IIPT,IGC) = EC(IIPT,IGC) + SCALE * DFDAB
C     -- DEBUG DENSITY
      RHOCHECK=RHOCHECK+RHOT*XYZWGT(IPT)
      E_LYP=E_LYP+ELYP*SCALE*XYZWGT(IPT)
      IF(NDER.EQ.1) GOTO 100
C
C     2ND DERIVATIVE
C
C      BASIC VALUABLES
C      ---------------
      TDDOMEGA = (C4/C3*OMEGA - RHOT*DOMEGA)*(C11*RHO13 - ALYPC
     * - ALYPD/XDNM)- RHOT*OMEGA* (C11/C3 - C1/C3*RX*RX)/RHO13/RHO13
      DDOMEGA = 1.0D+00/C3/RHOT/RHO43* TDDOMEGA
      TDDDELTA = -5.0D+00/C3/RHO43/RHO43/XDNM/XDNM +
     * C2/C3*ALYPD/RHO2/RHOT/XDNM/XDNM/XDNM
      DDDELTA = 1.0D+00/C3*((ALYPD*ALYPD)*TDDDELTA + 1.0D+00/RHO2
     * *DELTA-1.0D+00/RHOT*DDELTA)
      DOMEGA2 = -C4/C3/RHOT*(DOMEGA1 + 1.0D+00/C4*(C11/C3/RHOT
     *    - 1.0D+00/C3*ALYPD*ALYPD/XDNM/XDNM/RHO43/RHO13))
C      D2EDRA2 PART
C      ------------
      D2RA21A = (-1.0D+00/C3/XDNM/XDNM)*
     *          (C4/C3*ALYPD/RHOT/RHO43 + ALYPD*ALYPD/RHO43/RHO43)
     *         - 1.0D+00/RHOA/RHOA + 1.0D+00/RHOT/RHOT
      D2RA21_A = (-1.0D+00/C3/XDNM/XDNM)*
     *           (C4/C3*ALYPD/RHOT/RHO43 + ALYPD*ALYPD/RHO43/RHO43)
     *         - 1.0D+00/RHOB/RHOB + 1.0D+00/RHOT/RHOT
      D2RA21 = TERM2 * D2RA21A
      D2RA26AA = 1.0D+00/C3*RX*RHOA+RHOB
      D2RA26A =-C4*ALYPA/XDNM*RHOB/RHOT/RHOT*D2RA26AA
      D2RA26_AA = 1.0D+00/C3*RX*RHOB+RHOA
      D2RA26_A =-C4*ALYPA/XDNM*RHOA/RHOT/RHOT*D2RA26_AA
      D2RA26 = MBTERM2A * D2RA26A
      D2RA22A = (DDOMEGA + DOMEGA/RHOA)*TERM3
      D2RA22B =  C8/C3*DOMEGA*RHO83A*RHOB
      D2RA22C =  RHOB*(DOMEGA*TERM6A + C11*C8/C9*OMEGA*RHO83A/RHOA)
      D2RA22 = -(D2RA22A+D2RA22B+D2RA22C) * FACE * ALYPAB
      D2RA23A = (C2*RHOB/RHO2*DDELTA
     *         +(C3 + RHOA/RHOT)*DDDELTA
     *         - C2*(DELTA - C11)*RHOB/RHO2/RHOT)
      D2RA23_A = (C2*RHOA/RHO2*DDELTA
     *                 +(C3 + RHOB/RHOT)*DDDELTA
     *                 - C2*(DELTA - C11)*RHOA/RHO2/RHOT)
      D2RA23 = (DOMEGA2 * DFDAA + DOMEGA1 * D2AFDAA
     *    - ALYPAB * DOMEGA * MBABOAAA
     *    - ABOMEGA * (1.0D+00/C9*RHOB*(-C2*(TERM4A + TERM5B)
     *    -RHOA*D2RA23A))) * GRDAA
      D2RA24 = GRDAB*(DOMEGA2 * DFDAB + DOMEGA1 * D2AFDAB - ALYPAB
     * * DOMEGA * MBABOAAB- ABOMEGA * (-C7/C9*RHOB*(C2*DDELTA +
     * RHOA*DDDELTA) - C8/C3))
      D2RA25AA = (C3 + RHOB/RHOT)*DDDELTA
     *          -C2*RHOB/RHO2*DDELTA
     *          +C2*(DELTA-C11)*RHOB/RHOT/RHO2
      D2RA25A = -RHOB/C9*(C2*(TERM4B - TERM5B)
     *                    +RHOA*D2RA25AA)-C2
      D2RA25_AA = (C3 + RHOA/RHOT)*DDDELTA
     *                           -C2*RHOA/RHO2*DDELTA
     *                           +C2*(DELTA-C11)*RHOA/RHOT/RHO2
      D2RA25_A = -RHOA/C9*(C2*(TERM4A - TERM5A)
     *                    +RHOB*D2RA25_AA)-C2
      D2RA25 =  GRDBB*(DOMEGA2 * DFDBB + DOMEGA1 * D2AFDBB
     *         - ALYPAB * DOMEGA * MBABOABB
     *    - ABOMEGA * D2RA25A)
      D2EDRA2= D2RA21 + D2RA22 + D2RA23 + D2RA24 + D2RA25 +D2RA26
C      D2EDRB2 PART
C      ------------
      D2RB21 = TERM2 * ((-1.0D+00/C3/XDNM/XDNM)*(C4/C3*ALYPD/RHOT/RHO43
     *    + ALYPD*ALYPD/RHO43/RHO43) - 1.0D+00/RHOB/RHOB + 1.0D+00
     *    /RHOT/RHOT)
      D2RB26 =  MBTERM2B * (-C4*ALYPA/XDNM*RHOA/RHOT/RHOT*
     * (1.0D+00/C3*RX*RHOB+RHOA))
      D2RB22 = (DDOMEGA + DOMEGA/RHOB)*TERM3 + C8/C3*DOMEGA*RHO83B*RHOA
     *    + RHOA*(DOMEGA*TERM6B + C11*C8/C9*OMEGA*RHO83B/RHOB)
      D2RB22 = -D2RB22 * FACE * ALYPAB
      D2RB23 = (DOMEGA2 * DFDBB + DOMEGA1 * D2BFDBB
     *    - ABOMEGA * (1.0D+00/C9*RHOA*(-C2*(TERM4B + TERM5A)
     *    - RHOB*(C2*RHOA/RHO2*DDELTA + (C3 + RHOB/RHOT)*DDDELTA
     *    - C2*(DELTA - C11)*RHOA/RHO2/RHOT)))
     *    - ALYPAB * DOMEGA * MBABOBBB) * GRDBB
      D2RB24 = (DOMEGA2 * DFDAB + DOMEGA1 * D2BFDAB
     *    - ABOMEGA * (-C7/C9*RHOA*(C2*DDELTA + RHOB*DDDELTA) - C8/C3)
     *    - ALYPAB * DOMEGA * MBABOBAB) * GRDAB
      D2RB25 = (DOMEGA2 * DFDAA + DOMEGA1 * D2BFDAA
     *    - ABOMEGA * (-RHOA/C9*(C2*(TERM4A - TERM5A)
     *    + RHOB*((C3 + RHOA/RHOT)*DDDELTA - C2*RHOA/RHO2*DDELTA
     *    + C2*(DELTA-C11)*RHOA/RHOT/RHO2))-C2)
     *    - ALYPAB * DOMEGA * MBABOBAA) * GRDAA
      D2EDRB2= D2RB21 + D2RB22 + D2RB23 + D2RB24 + D2RB25 + D2RB26
C      D2EDRAB PART
C      ------------
      D2RAB1 = TERM2 * (-1.0D+00/C3/XDNM/XDNM*(C4/C3*ALYPD/RHOT/RHO43
     *    + ALYPD*ALYPD/RHO43/RHO43) + 1.0D+00/RHO2)
     *    + MBTERM2A * (-C4*ALYPA/XDNM*RHOA/RHO2*(1.0D+00/C3*RX*
     *    RHOB+RHOA))
      D2RAB2 = (DDOMEGA + DOMEGA/RHOB)*TERM3 + C8/C3*DOMEGA*RHOA*RHO83B
     *    + (DOMEGA*RHOB + OMEGA)*TERM6A + C8/C3*OMEGA*RHO83B
      D2RAB2 = -D2RAB2 * FACE * ALYPAB
      D2RAB3 = (DOMEGA2 * DFDAA + DOMEGA1 * D2BFDAA
     *    - ABOMEGA * (TERM1A + 1.0D+00/C9*RHOB*(-TERM4A + TERM5A)
     *    - 1.0D+00/C9*RHOA*(TERM4A + TERM5B)
     *    - 1.0D+00/C9*RHOA*RHOB*((RHOB-RHOA)/RHO2*DDELTA +
     *    (C3 + RHOA/RHOT)*DDDELTA
     *    + (DELTA-C11)*(RHO2A - RHO2B)/RHO2/RHO2))
     *    - ALYPAB * DOMEGA * MBABOAAA) * GRDAA
      D2RAB4 = (DOMEGA2 * DFDAB + DOMEGA1 * D2BFDAB
     *    - ABOMEGA * (1.0D+00/C9*(C47-C7*DELTA)-C7/C9*(RHOT*DDELTA
     *    + RHOA*RHOB*DDDELTA)
     *    - C8/C3) - ALYPAB * DOMEGA * MBABOAAB) * GRDAB
      D2RAB5 = (DOMEGA2 * DFDBB + DOMEGA1 * D2BFDBB
     *    - ABOMEGA * (TERM1B + 1.0D+00/C9*RHOB*(-TERM4B - TERM5A)
     *    - 1.0D+00/C9*RHOA*(TERM4B - TERM5B)
     *    - 1.0D+00/C9*RHOA*RHOB*((C3 + RHOB/RHOT)*DDDELTA
     *    + (RHOA - RHOB)/RHO2*DDELTA
     *    - (DELTA-C11)*(RHO2A - RHO2B)/RHO2/RHO2))
     *    - ALYPAB * DOMEGA * MBABOABB) * GRDBB
      D2EDRAB= D2RAB1 + D2RAB2 + D2RAB3 + D2RAB4 + D2RAB5
C     2ND DERIVATIVE
      EC(IIPT,IRARA)=EC(IIPT,IRARA) + SCALE*D2EDRA2
      EC(IIPT,IRBRB)=EC(IIPT,IRBRB) + SCALE*D2EDRB2
      EC(IIPT,IRARB)=EC(IIPT,IRARB) + SCALE*D2EDRAB
      EC(IIPT,IRAGA)=EC(IIPT,IRAGA) + SCALE*D2AFDAA
      EC(IIPT,IRBGA)=EC(IIPT,IRBGA) + SCALE*D2BFDAA
      EC(IIPT,IRAGB)=EC(IIPT,IRAGB) + SCALE*D2AFDBB
      EC(IIPT,IRBGB)=EC(IIPT,IRBGB) + SCALE*D2BFDBB
      EC(IIPT,IRAGC)=EC(IIPT,IRAGC) + SCALE*D2AFDAB
      EC(IIPT,IRBGC)=EC(IIPT,IRBGC) + SCALE*D2BFDAB
      IF(NDER.EQ.2) GOTO 100
C
C      THIRD DERIVATIVES
C     --------------------------------
C     DXC    =-ALYPC/(C3*RHO43)
      DXD    =-ALYPD/(C3*RHO43)
      DXDNM  = DXD
      DRX    = -ALYPD/(C3*(ALYPD**2*RHO23+RHO43+C2*ALYPD*RHOT))
      ABDOMEGA=ALYPAB*DOMEGA
      DOMEGA3 = C4/C3/RHO2*(-RHOT*DOMEGA2+DOMEGA1+C11/6/RHOT
     * +ALYPD*ALYPD/18/XDNM/XDNM*(ALYPD/XDNM/RHO2-4/RHO53))
      D3OMEGA = -C7/C9/RHO53/RHO53*TDDOMEGA + C1/C3/RHO73*(
     * (RHOT*DDOMEGA-DOMEGA/C3)*(ALYPC+ALYPD/XDNM-C11*RHO13)
     * +(C11/C3/RHO23+ALYPD*DXDNM/XDNM/XDNM)*(C4/C3*OMEGA-DOMEGA*RHOT)
     * +((RX*RX-C11)*(OMEGA+RHOT*DOMEGA)/C3+C2/C3*RHOT*RX*DRX*OMEGA)
     * /RHO23 -C2/C3/RHO53*(RHOT*OMEGA*(-C11/C3+RX*RX/C3)))
      D3DELTA = (ALYPD*ALYPD*(C40/C9/RHO113/XDNM/XDNM-C28*ALYPD/C9/
     * RHOT**4/XDNM**3+C2*ALYPD*ALYPD/C3/RHO133/XDNM**4)
     *          +C2*DDELTA/RHO2-DDDELTA/RHOT-C2*DELTA/RHO2/RHOT) /C3
      DTERM1ADA =-(DDELTA*(3*RHO2+RHOA*RHOT)+RHOB*DELTA-11*RHOB)/C9/RHO2
      DTERM1ADB =-(DDELTA*(3*RHO2+RHOA*RHOT)-RHOA*DELTA+11*RHOA)/C9/RHO2
      DTERM1BDA =-(DDELTA*(3*RHO2+RHOB*RHOT)-RHOB*DELTA+11*RHOB)/C9/RHO2
      DTERM1BDB =-(DDELTA*(3*RHO2+RHOB*RHOT)+RHOA*DELTA-11*RHOA)/C9/RHO2
      DTERM2DA = -C4*ALYPA*RHOB*(-C3*(RHOA*RHO13-RHO43)+ALYPD*
     * (RHOT+2*RHOB))/(C3*(ALYPD*ALYPD*RHO53+RHO73+2*ALYPD*RHO2))
      DTERM2DB = -C4*ALYPA*RHOA*(-C3*(RHOB*RHO13-RHO43)+ALYPD*
     * (RHOT+2*RHOA))/(C3*(ALYPD*ALYPD*RHO53+RHO73+2*ALYPD*RHO2))
      DTERM3DA = RHOB*(C11/C3*RHO83A+RHO83B)
      DTERM3DB = RHOA*(C11/C3*RHO83B+RHO83A)
      DTERM4ADA = (DDDELTA*(C4*RHO2A+C3*RHO2B+C7*RHOA*RHOB)+RHOB*
     * DDELTA)/RHO2
      DTERM4ADB = (DDDELTA*(C4*RHO2A+C3*RHO2B+C7*RHOA*RHOB)-RHOA*
     * DDELTA)/RHO2
      DTERM4BDA = (DDDELTA*(C4*RHO2B+C3*RHO2A+C7*RHOA*RHOB)-RHOB*
     * DDELTA)/RHO2
      DTERM4BDB = (DDDELTA*(C4*RHO2B+C3*RHO2A+C7*RHOA*RHOB)+RHOA*
     * DDELTA)/RHO2
      DTERM5ADA = ((DELTA-C11)*(RHOB-RHOA)+RHOA*DDELTA*RHOT)/RHO2/RHOT
      DTERM5ADB = RHOA*(-2*DELTA+RHOT*DDELTA+22)/RHO2/RHOT
      DTERM5BDA = RHOB*(-2*DELTA+RHOT*DDELTA+22)/RHO2/RHOT
      DTERM5BDB = ((DELTA-C11)*(RHOA-RHOB)+RHOB*DDELTA*RHOT)/RHO2/RHOT
      DTERM6ADA = 88*RHO43A*RHO13A/C9
      DTERM6ADB = C8*RHO43B*RHO13B/C3
      DTERM6BDA = C8*RHO43A*RHO13A/C3
      DTERM6BDB = 88*RHO43B*RHO13B/C9
      MBABOAAADA = RHOB*DTERM1ADA-RHOB/C9*(TERM4A+TERM5B)-RHOA*RHOB/C9
     * *(DTERM4ADA+DTERM5BDA)
      MBABOAAADB = TERM1A+RHOB*DTERM1ADB-RHOA/C9*(TERM4A+TERM5B)-RHOA
     * *RHOB/C9*(DTERM4ADB+DTERM5BDB)
      MBABOABBDA = RHOB*DTERM1BDA -RHOB/C9*(TERM4B-TERM5B)-RHOA*RHOB/
     * C9*(DTERM4BDA-DTERM5BDA)-C2
      MBABOABBDB = TERM1B+RHOB*DTERM1BDB-RHOA/C9*(TERM4B-TERM5B)
     * -RHOA*RHOB/C9*(DTERM4BDB-DTERM5BDB)
      MBABOAABDA = -RHOB/C9*C7*DDELTA - C7/C9*RHOB*(RHOA*DDDELTA+DDELTA)
     * -C8/C3
      MBABOAABDB = ((C47-C7*DELTA)-C7*RHOB*DDELTA)/C9-C7/C9*(RHOA*
     * DDELTA+RHOA*RHOB*DDDELTA)- C8/C3
      MBABOBBBDA = TERM1B+RHOA*DTERM1BDA-RHOB/C9*(TERM4B+TERM5A)
     * -RHOA*RHOB/C9*(DTERM4BDA+DTERM5ADA)
      MBABOBBBDB = RHOA*DTERM1BDB-RHOA/C9*(TERM4B+TERM5A)-RHOA*
     * RHOB/C9*(DTERM4BDB+DTERM5ADB)
      MBABOBAADA = TERM1A+RHOA*DTERM1ADA -RHOB/C9*(TERM4A-TERM5A)
     * -RHOA*RHOB/C9*(DTERM4ADA-DTERM5ADA)
      MBABOBAADB = RHOA*DTERM1ADB        -RHOA/C9*(TERM4A-TERM5A)
     * -RHOA*RHOB/C9*(DTERM4ADB-DTERM5ADB)-C2
      MBABOBABDA = (C47-C7*DELTA)/C9+RHOA*(-C7*DDELTA)/C9-C7/C9*
     * RHOB*(DDELTA+RHOA*DDDELTA)-C8/C3
      MBABOBABDB = RHOA*(-C7*DDELTA)/C9-C7/C9*RHOA*(DDELTA+RHOB*
     * DDDELTA)-C8/C3
      MBTERM2ADB =-ALYPD*C4/C9/RHO73/XDNM-ALYPD/RHO43/C3*DXDNM/XDNM/XDNM
     *           +C1/RHO2
      MBTERM2BDA = MBTERM2ADB
      MBTERM2ADA = MBTERM2ADB-C1/RHO2A
      MBTERM2BDB = MBTERM2ADB-C1/RHO2B
      D3EDRARAGA =
     * DOMEGA/OMEGA*D2AFDAA+DFDAA*(DDOMEGA/OMEGA-(DOMEGA/OMEGA)**2)
     * - ALYPAB*(
     * DOMEGA*MBABOAAA
     * -RHOB*OMEGA/C9*(2*((DELTA-C11)*RHOB/RHO2+(C3+RHOA/RHOT)*DDELTA)
     * +RHOA*(2*RHOB/RHO2*DDELTA+(C3+RHOA/RHOT)*DDDELTA
     * -(DELTA-C11)*(C2*RHOB)/(RHOT**3))))
      D3EDRBRBGB =
     * DOMEGA/OMEGA*D2BFDBB+DFDBB*(DDOMEGA/OMEGA-(DOMEGA/OMEGA)**2)
     * - ALYPAB*(
     * DOMEGA*MBABOBBB
     * -RHOA*OMEGA/C9*(2*((DELTA-C11)*RHOA/RHO2+(C3+RHOB/RHOT)*DDELTA)
     * +RHOB*(2*RHOA/RHO2*DDELTA+(C3+RHOB/RHOT)*DDDELTA
     * -(DELTA-C11)*(C2*RHOA)/(RHOT**3))))
      D3EDRARBGA =
     * DOMEGA/OMEGA*D2BFDAA+DFDAA*(DDOMEGA/OMEGA-(DOMEGA/OMEGA)**2)
     * - ALYPAB*(DOMEGA*MBABOAAA +OMEGA*MBABOAAADB)
      D3EDRARBGB =
     * DOMEGA/OMEGA*D2AFDBB+DFDBB*(DDOMEGA/OMEGA-(DOMEGA/OMEGA)**2)
     * - ALYPAB*(DOMEGA*MBABOBBB+OMEGA*MBABOBBBDA)
      D3EDRARAGC =
     * DOMEGA/OMEGA*D2AFDAB+DFDAB*(DDOMEGA/OMEGA-(DOMEGA/OMEGA)**2)
     * - ALYPAB*(DOMEGA*MBABOAAB +
     * OMEGA*(-C7/C9*RHOA*RHOB*DDDELTA-C8/C3-C2*C7/C9*RHOB*DDELTA))
      D3EDRBRBGC =
     * DOMEGA/OMEGA*D2BFDAB+DFDAB*(DDOMEGA/OMEGA-(DOMEGA/OMEGA)**2)
     * - ALYPAB*(DOMEGA*MBABOBAB +
     * OMEGA*(-C7/C9*RHOA*RHOB*DDDELTA-C8/C3-C2*C7/C9*RHOA*DDELTA))
      D3EDRARBGC =
     * DOMEGA/OMEGA*D2BFDAB+DFDAB*(DDOMEGA/OMEGA-(DOMEGA/OMEGA)**2)
     * - ALYPAB*(DOMEGA*MBABOAAB +
     * OMEGA*((47-C7*DELTA)/C9-C7/C9*RHOB*DDELTA-C7/C9*(RHOA*DDELTA
     * +RHOA*RHOB*DDDELTA)-C8/C3))
      D3EDRARAGB =
     * DOMEGA/OMEGA*D2AFDBB+DFDBB*(DDOMEGA/OMEGA-(DOMEGA/OMEGA)**2)
     * - ALYPAB*(DOMEGA*MBABOABB +
     * OMEGA*(RHOB/C9*(-C3*DDELTA-DDELTA*RHOB/RHOT+(DELTA-11)*RHOB/RHO2)
     * -RHOB/C9*((C3+RHOB/RHOT)*DDELTA-(DELTA-11)*RHOB/RHO2)
     * -RHOA*RHOB/C9*(-RHOB/RHO2*DDELTA+(3+RHOB/RHOT)*DDDELTA-DDELTA
     * *RHOB/RHO2+(DELTA-11)*2*RHOB/RHOT**3)-2))
      D3EDRBRBGA =
     * DOMEGA/OMEGA*D2BFDAA+DFDAA*(DDOMEGA/OMEGA-(DOMEGA/OMEGA)**2)
     * - ALYPAB*(DOMEGA*MBABOBAA +
     * OMEGA*(RHOA/C9*(-C3*DDELTA-DDELTA*RHOA/RHOT+(DELTA-11)*RHOA/RHO2)
     * -RHOA/C9*((C3+RHOA/RHOT)*DDELTA-(DELTA-11)*RHOA/RHO2)
     * -RHOB*RHOA/C9*(-RHOA/RHO2*DDELTA+(3+RHOA/RHOT)*DDDELTA
     * -DDELTA*RHOA/RHO2+(DELTA-11)*2*RHOA/RHOT**3)-2))
      D3RA21A = TERM2*((C8*ALYPD2/C3/RHO113+C28*ALYPD/C9/RHO103)
     * /C3/XDNM/XDNM
     * +C2*DXDNM*(ALYPD2/RHO83+C4*ALYPD/C3/RHO73)/C3/XDNM/XDNM/XDNM
     * -C2/RHO2/RHOT+C2/RHO2A/RHOA) +DTERM2DA*D2RA21A
      D3RA21B = TERM2*((C8*ALYPD2/C3/RHO113+C28*ALYPD/C9/RHO103)/C3
     * /XDNM/XDNM
     * +C2*DXDNM*(ALYPD2/RHO83+C4*ALYPD/C3/RHO73)/C3/XDNM/XDNM/XDNM
     * -C2/RHO2/RHOT) +DTERM2DB*D2RA21A
      D3RB21B = TERM2*((C8*ALYPD2/C3/RHO113+C28*ALYPD/C9/RHO103)
     * /C3/XDNM/XDNM
     * +C2*DXDNM*(ALYPD2/RHO83+C4*ALYPD/C3/RHO73)/C3/XDNM/XDNM/XDNM
     * -C2/RHO2/RHOT+C2/RHO2B/RHOB) +DTERM2DB*D2RA21_A
      D3RB21A = TERM2*((C8*ALYPD2/C3/RHO113+C28*ALYPD/C9/RHO103)
     * /C3/XDNM/XDNM
     * +C2*DXDNM*(ALYPD2/RHO83+C4*ALYPD/C3/RHO73)/C3/XDNM/XDNM/XDNM
     * -C2/RHO2/RHOT) +DTERM2DA*D2RA21_A
      D3RA22A = TERM3*(D3OMEGA-DOMEGA/RHO2A+DDOMEGA/RHOA)+DTERM3DA
     * *(DDOMEGA+DOMEGA/RHOA)
     * +C8/C3*RHOB*(DDOMEGA*RHO83A+C8/C3*DOMEGA*RHO43A*RHO13A)
     * +RHOB*(DDOMEGA*TERM6A+DOMEGA*DTERM6ADA
     * + C11*C8/C9*(DOMEGA*RHO83A/RHOA+C5/C3*OMEGA*RHO13A*RHO13A))
      D3RA22A = - D3RA22A * FACE*ALYPAB
      D3RA22B = TERM3*(D3OMEGA+DDOMEGA/RHOA)+DTERM3DB*
     * (DDOMEGA+DOMEGA/RHOA)
     * +C8/C3*RHO83A*(DOMEGA+RHOB*DDOMEGA)
     * +(DOMEGA*TERM6A + C11*C8/C9*OMEGA*RHO83A/RHOA)
     * +RHOB*(DDOMEGA*TERM6A+DOMEGA*DTERM6ADB
     * +C11*C8/C9*RHO83A/RHOA*DOMEGA)
      D3RA22B = - D3RA22B * FACE*ALYPAB
      D3RB22B = TERM3*(D3OMEGA-DOMEGA/RHO2B+DDOMEGA/RHOB)
     * +DTERM3DB*(DDOMEGA+DOMEGA/RHOB)
     * +C8/C3*RHOA*(DDOMEGA*RHO83B+C8/C3*DOMEGA*RHO43B*RHO13B)
     * +RHOA*(DDOMEGA*TERM6B+DOMEGA*DTERM6BDB
     * + C11*C8/C9*(DOMEGA*RHO83B/RHOB+C5/C3*OMEGA*RHO13B*RHO13B))
      D3RB22B = - D3RB22B * FACE*ALYPAB
      D3RB22A = TERM3*(D3OMEGA+DDOMEGA/RHOB)+DTERM3DA*
     * (DDOMEGA+DOMEGA/RHOB)
     * +C8/C3*RHO83B*(DOMEGA+RHOA*DDOMEGA)
     * +(DOMEGA*TERM6B + C11*C8/C9*OMEGA*RHO83B/RHOB)
     * +RHOA*(DDOMEGA*TERM6B+DOMEGA*DTERM6BDA
     * +C11*C8/C9*RHO83B/RHOB*DOMEGA)
      D3RB22A = - D3RB22A * FACE*ALYPAB
      D3RA23A = GRDAA*(DOMEGA3*DFDAA+DOMEGA2*D2AFDAA+DOMEGA2*
     * D2AFDAA+DOMEGA1*D3EDRARAGA
     * -ALYPAB*(DDOMEGA*MBABOAAA+DOMEGA*MBABOAAADA)
     * -ABDOMEGA*(C1/C9*RHOB*(-C2*(TERM4A+TERM5B)-RHOA*D2RA23A))
     * -ABOMEGA*(RHOB/C9*(-C2*(DTERM4ADA+DTERM5BDA)-D2RA23A
     * -RHOA*(C2*RHOB*(-C2/RHO3*DDELTA+DDDELTA/RHO2)
     * +(C3+RHOA/RHOT)*D3DELTA+(1/RHOT-RHOA/RHO2)*DDDELTA
     * -C2*RHOB*(DDELTA/RHO3-C3*(DELTA-C11)/RHO2/RHO2)))))
      D3RB23B = GRDBB*(DOMEGA3*DFDBB+DOMEGA2*D2BFDBB+DOMEGA2*D2BFDBB
     * +DOMEGA1*D3EDRBRBGB
     * -ALYPAB*(DDOMEGA*MBABOBBB+DOMEGA*MBABOBBBDB)
     * -ABDOMEGA*(C1/C9*RHOA*(-C2*(TERM4B+TERM5A)-RHOB*D2RA23_A))
     * -ABOMEGA*(RHOA/C9*(-C2*(DTERM4BDB+DTERM5ADB)-D2RA23_A
     * -RHOB*(C2*RHOA*(-C2/RHO3*DDELTA+DDDELTA/RHO2)
     * +(C3+RHOB/RHOT)*D3DELTA+(1/RHOT-RHOB/RHO2)*DDDELTA
     * -C2*RHOA*(DDELTA/RHO3-C3*(DELTA-C11)/RHO2/RHO2)))))
      D2RA23B = (C2/RHO2-C4*RHOB/RHO3)*DDELTA+(C2*RHOB/RHO2*DDDELTA)
     * +(-RHOA/RHO2)*DDDELTA + (C3+RHOA/RHOT)*D3DELTA
     * -C2*(DDELTA*RHOB/RHO3 + (DELTA-C11)*(C1/RHO3-C3*RHOB/RHO2/RHO2))
      D3RA23B = GRDAA*(DOMEGA3*DFDAA+DOMEGA2*D2BFDAA+DOMEGA2*D2AFDAA+
     * DOMEGA1*D3EDRARBGA
     * -ALYPAB*(DDOMEGA*MBABOAAA+DOMEGA*MBABOAAADB)
     * -ABDOMEGA*(C1/C9*RHOB*(-C2*(TERM4A+TERM5B)-RHOA*D2RA23A))
     * -ABOMEGA/C9*(-C2*(TERM4A + TERM5B)-RHOA*D2RA23A
     * +RHOB*(-C2*(DTERM4ADB+DTERM5BDB)-RHOA*D2RA23B)))
      D3RB23A = GRDBB*(DOMEGA3*DFDBB+DOMEGA2*D2AFDBB+DOMEGA2*D2BFDBB
     * +DOMEGA1*D3EDRARBGB
     * -ALYPAB*(DDOMEGA*MBABOBBB+DOMEGA*MBABOBBBDA)
     * -ABDOMEGA*(C1/C9*RHOA*(-C2*(TERM4B+TERM5A)-RHOB*D2RA23_A))
     * -ABOMEGA/C9*(-C2*(TERM4B + TERM5A)-RHOB*D2RA23_A
     * +RHOA*(-C2*(DTERM4BDA+DTERM5ADA)-RHOB*(
     * (C2/RHO2-C4*RHOA/RHO3)*DDELTA+(C2*RHOA/RHO2*DDDELTA)
     * +(-RHOB/RHO2)*DDDELTA + (C3+RHOB/RHOT)*D3DELTA
     * -C2*(DDELTA*RHOA/RHO3 + (DELTA-C11)*(C1/RHO3-C3*RHOA/RHO2/RHO2))
     * ))))
      D3RA24A = GRDAB*(DOMEGA3*DFDAB+DOMEGA2*D2AFDAB
     * +DOMEGA2*D2AFDAB+DOMEGA1*D3EDRARAGC
     * -ALYPAB*DDOMEGA*MBABOAAB - ABDOMEGA*MBABOAABDA
     * -ABDOMEGA*(-C7/C9*RHOB*(C2*DDELTA+RHOA*DDDELTA)-C8/C3)
     * +ABOMEGA*C7/C9*RHOB*(C2*DDDELTA + DDDELTA + RHOA*D3DELTA))
      D3RB24B = GRDAB*(DOMEGA3*DFDAB+DOMEGA2*D2BFDAB
     * +DOMEGA2*D2BFDAB+DOMEGA1*D3EDRBRBGC
     * -ALYPAB*DDOMEGA*MBABOBAB - ABDOMEGA*MBABOBABDB
     * -ABDOMEGA*(-C7/C9*RHOA*(C2*DDELTA+RHOB*DDDELTA)-C8/C3)
     * +ABOMEGA*C7/C9*RHOA*(C2*DDDELTA + DDDELTA + RHOB*D3DELTA))
      D3RA24B = GRDAB*(DOMEGA3*DFDAB+DOMEGA2*D2BFDAB
     * +DOMEGA2*D2AFDAB+DOMEGA1*D3EDRARBGC
     * -ALYPAB*DDOMEGA*MBABOAAB - ABDOMEGA*MBABOAABDB
     * -ABDOMEGA*(-C7/C9*RHOB*(C2*DDELTA + RHOA*DDDELTA) - C8/C3)
     * -ABOMEGA*(-C7/C9*((C2*DDELTA + RHOA*DDDELTA)+RHOB*
     * (C2*DDDELTA+RHOA*D3DELTA))))
      D3RB24A = GRDAB*(DOMEGA3*DFDAB+DOMEGA2*D2AFDAB
     * +DOMEGA2*D2BFDAB+DOMEGA1*D3EDRARBGC
     * -ALYPAB*DDOMEGA*MBABOBAB - ABDOMEGA*MBABOBABDA
     * -ABDOMEGA*(-C7/C9*RHOA*(C2*DDELTA + RHOB*DDDELTA) - C8/C3)
     * -ABOMEGA*(-C7/C9*((C2*DDELTA + RHOB*DDDELTA)
     * +RHOA*(C2*DDDELTA+RHOB*D3DELTA))))
      D3RA25A = GRDBB*(DOMEGA3*DFDBB+DOMEGA2*D2AFDBB
     * +DOMEGA2*D2AFDBB+DOMEGA1*D3EDRARAGB
     * -ALYPAB*DDOMEGA*MBABOABB-ABDOMEGA*MBABOABBDA
     * -ABDOMEGA*D2RA25A - ABOMEGA*(
     * -RHOB/C9*(C2*(DTERM4BDA - DTERM5BDA)
     * +D2RA25AA+RHOA*(
     * -RHOB/RHO2*DDDELTA+(C3+RHOB/RHOT)*D3DELTA
     * +C4*RHOB/RHO3*DDELTA-C2*RHOB/RHO2*DDDELTA
     * +C2*DDELTA*RHOB/RHO3-C6*(DELTA-C11)*RHOB/RHO2/RHO2
     * ))))
      D3RB25B = GRDAA*(DOMEGA3*DFDAA+DOMEGA2*D2BFDAA
     * +DOMEGA2*D2BFDAA+DOMEGA1*D3EDRBRBGA
     * -ALYPAB*DDOMEGA*MBABOBAA-ABDOMEGA*MBABOBAADB
     * -ABDOMEGA*D2RA25_A - ABOMEGA*(
     * -RHOA/C9*(C2*(DTERM4ADB - DTERM5ADB)
     * +D2RA25_AA+RHOB*(
     * -RHOA/RHO2*DDDELTA+(C3+RHOA/RHOT)*D3DELTA
     * +C4*RHOA/RHO3*DDELTA-C2*RHOA/RHO2*DDDELTA
     * +C2*DDELTA*RHOA/RHO3-C6*(DELTA-C11)*RHOA/RHO2/RHO2
     * ))))
      D3RA25B = GRDBB*(DOMEGA3*DFDBB+DOMEGA2*D2BFDBB+DOMEGA2*D2AFDBB
     * +DOMEGA1*D3EDRARBGB
     * -ALYPAB*DDOMEGA*MBABOABB-ABDOMEGA*MBABOABBDB
     * -ABDOMEGA*D2RA25A-ABOMEGA*(
     * -C1/C9*(C2*(TERM4B-TERM5B)+RHOA*D2RA25AA)
     * -RHOB/C9*(C2*(DTERM4BDB-DTERM5BDB)+RHOA*(
     * (C1/RHOT-RHOB/RHO2)*DDDELTA+(C3+RHOB/RHOT)*D3DELTA
     * -C2*(C1/RHO2-C2*RHOB/RHO3)*DDELTA-C2*RHOB/RHO2*DDDELTA
     * +C2*(DELTA-C11)*(C1/RHO3-C3*RHOB/RHO2/RHO2)+C2*DDELTA*RHOB/RHO3)
     * )))
      D3RB25A = GRDAA*(DOMEGA3*DFDAA+DOMEGA2*D2AFDAA+DOMEGA2*D2BFDAA
     * +DOMEGA1*D3EDRARBGA
     * -ALYPAB*DDOMEGA*MBABOBAA-ABDOMEGA*MBABOBAADA
     * -ABDOMEGA*D2RA25_A-ABOMEGA*(
     * -C1/C9*(C2*(TERM4A-TERM5A)+RHOB*D2RA25_AA)
     * -RHOA/C9*(C2*(DTERM4ADA-DTERM5ADA)+RHOB*(
     * (C1/RHOT-RHOA/RHO2)*DDDELTA+(C3+RHOA/RHOT)*D3DELTA
     * -C2*(C1/RHO2-C2*RHOA/RHO3)*DDELTA-C2*RHOA/RHO2*DDDELTA
     * +C2*(DELTA-C11)*(C1/RHO3-C3*RHOA/RHO2/RHO2)+C2*DDELTA*RHOA/RHO3)
     * )))
      D3RA26A =  MBTERM2ADA*D2RA26A+MBTERM2A*(
     *            -C4*ALYPA*RHOB*( -DXDNM/XDNM2/RHO2*D2RA26AA
     *                           -C2*D2RA26AA/RHO3/XDNM
     *                           +(DRX*RHOA+RX)/RHO2/C3/XDNM))
      D3RB26B =  MBTERM2BDB*D2RA26_A+MBTERM2B*(
     *            -C4*ALYPA*RHOA*( -DXDNM/XDNM2/RHO2*D2RA26_AA
     *                           -C2*D2RA26_AA/RHO3/XDNM
     *                           +(DRX*RHOB+RX)/RHO2/C3/XDNM))
      D3RA26B =  MBTERM2ADB*D2RA26A+MBTERM2A*(
     *            -C4*ALYPA*( -DXDNM/XDNM2/RHO2*D2RA26AA*RHOB
     *                       -C2*D2RA26AA/RHO3/XDNM*RHOB
     *                       +(RHOA*DRX/C3+C1)/RHO2/XDNM*RHOB
     *                       +D2RA26AA/RHO2/XDNM))
      D3RB26A =  MBTERM2BDA*D2RA26_A+MBTERM2B*(
     *            -C4*ALYPA*( -DXDNM/XDNM2/RHO2*D2RA26_AA*RHOA
     *                       -C2*D2RA26_AA/RHO3/XDNM*RHOA
     *                       +(RHOB*DRX/C3+C1)/RHO2/XDNM*RHOA
     *                       +D2RA26_AA/RHO2/XDNM))
C
      D3EDRARARA =  D3RA21A + D3RA22A + D3RA23A + D3RA24A
     * + D3RA25A +D3RA26A
      D3EDRARARB =  D3RA21B + D3RA22B + D3RA23B + D3RA24B + D3RA25B
     * +D3RA26B
      D3EDRBRBRB =  D3RB21B + D3RB22B + D3RB23B + D3RB24B + D3RB25B
     * +D3RB26B
      D3EDRARBRB =  D3RB21A + D3RB22A + D3RB23A + D3RB24A + D3RB25A
     * +D3RB26A
C
      EC(IIPT,IRARARA) = EC(IIPT,IRARARA) + SCALE*D3EDRARARA
      EC(IIPT,IRARARB) = EC(IIPT,IRARARB) + SCALE*D3EDRARARB
      EC(IIPT,IRARAGA) = EC(IIPT,IRARAGA) + SCALE*D3EDRARAGA
      EC(IIPT,IRARAGB) = EC(IIPT,IRARAGB) + SCALE*D3EDRARAGB
      EC(IIPT,IRARAGC) = EC(IIPT,IRARAGC) + SCALE*D3EDRARAGC
      EC(IIPT,IRARBRB) = EC(IIPT,IRARBRB) + SCALE*D3EDRARBRB
      EC(IIPT,IRARBGA) = EC(IIPT,IRARBGA) + SCALE*D3EDRARBGA
      EC(IIPT,IRARBGB) = EC(IIPT,IRARBGB) + SCALE*D3EDRARBGB
      EC(IIPT,IRARBGC) = EC(IIPT,IRARBGC) + SCALE*D3EDRARBGC
      EC(IIPT,IRBRBRB) = EC(IIPT,IRBRBRB) + SCALE*D3EDRBRBRB
      EC(IIPT,IRBRBGA) = EC(IIPT,IRBRBGA) + SCALE*D3EDRBRBGA
      EC(IIPT,IRBRBGB) = EC(IIPT,IRBRBGB) + SCALE*D3EDRBRBGB
      EC(IIPT,IRBRBGC) = EC(IIPT,IRBRBGC) + SCALE*D3EDRBRBGC
      ENDIF
  100 CONTINUE
      END DO
      ECOR=ECOR+E_LYP
C
      IF(DBUG) THEN
      WRITE(6,*)'INTEGRATED RHO IN DERLYP=',RHOCHECK
      WRITE(6,*)'LYP CORRELATION ENERGY=',E_LYP
      ENDIF
C
      RETURN
      END
C
C     INTERFACING THE EXISTING DFT FUNCTIONALS WHICH CONTAINS 2ND
C     FUNCTIONAL DERIVATIVES WITH THE TDDFT ENERGY CODE
C
C*MODULE TDFFUN  *DECK NDFPW91CF
      SUBROUTINE NDFPW91CF(ECOR,EC0,EC,DRHOI_,XYZWGT,SCALE,RHOMIN,
     *                     ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      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 EC(ILENG,*)
      DIMENSION EC0(*)
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION XYZWGT(*)
      DATA ZERO/0.0D+00/
      DATA TWO /2.0D+00/
C
      IDERIV = NDER
C     IDERIV = NUMBER OF DERIVATIVES NEEDED
C              TDDFT ENERGIES REQUIRE NDER=2
C              TDDFT GRADIENTS REQUIRE NDER=3
C
      NPT    = 1
C     NPT    = NUMBER POINTS TO BE EVALUATED IN THE DFT REPOSITORY
C              ROUTINE
C
      ECF    = ZERO
C     ECF    = INITIALIZE THE ENERGY OF THE CORRELATION FUNCTION TO ZERO
C
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
         RHOA=DRHOI_(IPT,1,1)
         RHOB=DRHOI_(IPT,1,IRU)
C             IN THE CLOSED SHELL CASE RHOB =RHOA SINCE IRU=1
C             IN THE OPENED SHELL CASE RHOB!=RHOA SINCE IRU=2
         RHOT=RHOA + RHOB
C
         IF(RHOT/TWO.GT.RHOMIN) THEN
           DXA = DRHOI_(IPT,2,1)
           DYA = DRHOI_(IPT,3,1)
           DZA = DRHOI_(IPT,4,1)
           DXB = DRHOI_(IPT,2,IRU)
           DYB = DRHOI_(IPT,3,IRU)
           DZB = DRHOI_(IPT,4,IRU)
           GRDAA=DXA*DXA + DYA*DYA + DZA*DZA
           GRDBB=DXB*DXB + DYB*DYB + DZB*DZB
           GRDAB=DXA*DXB + DYA*DYB + DZA*DZB
C
           RHOA1=RHOA
           RHOB1=RHOB
           SIGMAAA1=GRDAA
           SIGMABB1=GRDBB
           SIGMAAB1=GRDAB
C
C          INITIALIZE FIRST DERIVATIVES
C
           VRHOA                = 0.00D+00
           VRHOB                = 0.00D+00
           VSIGMAAA             = 0.00D+00
           VSIGMABB             = 0.00D+00
           VSIGMAAB             = 0.00D+00
C
C          INITIALIZE SECOND DERIVATIVES
C
           V2RHOA2              = 0.00D+00
           V2RHOAB              = 0.00D+00
           V2RHOASIGMAAA        = 0.00D+00
           V2RHOASIGMABB        = 0.00D+00
           V2RHOASIGMAAB        = 0.00D+00
           V2RHOB2              = 0.00D+00
           V2RHOBSIGMAAA        = 0.00D+00
           V2RHOBSIGMABB        = 0.00D+00
           V2RHOBSIGMAAB        = 0.00D+00
           V2SIGMAAA2           = 0.00D+00
           V2SIGMAAABB          = 0.00D+00
           V2SIGMAAAAB          = 0.00D+00
           V2SIGMABB2           = 0.00D+00
           V2SIGMAABBB          = 0.00D+00
           V2SIGMAAB2           = 0.00D+00
C
           CALL UKS_C_PW91
     &    (IDERIV,NPT,RHOA1,RHOB1,SIGMAAA1,SIGMABB1,SIGMAAB1,
     &     ZK,VRHOA,VRHOB,VSIGMAAA,VSIGMABB,VSIGMAAB,
     &     V2RHOA2,V2RHOB2,V2RHOAB,
     &     V2RHOASIGMAAA,V2RHOASIGMAAB,V2RHOASIGMABB,
     &     V2RHOBSIGMABB,V2RHOBSIGMAAB,V2RHOBSIGMAAA,
     &     V2SIGMAAA2,V2SIGMAAAAB,V2SIGMAAABB,
     &     V2SIGMAAB2,V2SIGMAABBB,V2SIGMABB2)
C
            ECF             = ECF             + SCALE*ZK*XYZWGT(IPT)
            EC0(IIPT)       = EC0(IIPT)       + SCALE*ZK
C           FIRST DERIVATIVES
            EC (IIPT,IRA)   = EC (IIPT,IRA)   + SCALE*VRHOA
            EC (IIPT,IRB)   = EC (IIPT,IRB)   + SCALE*VRHOB
            EC (IIPT,IGA)   = EC (IIPT,IGA)   + SCALE*VSIGMAAA
            EC (IIPT,IGB)   = EC (IIPT,IGB)   + SCALE*VSIGMABB
            EC (IIPT,IGC)   = EC (IIPT,IGC)   + SCALE*VSIGMAAB
C           SECOND DERIVATIVES
            EC (IIPT,IRARA) = EC (IIPT,IRARA) + SCALE*V2RHOA2
            EC (IIPT,IRARB) = EC (IIPT,IRARB) + SCALE*V2RHOAB
            EC (IIPT,IRAGA) = EC (IIPT,IRAGA) + SCALE*V2RHOASIGMAAA
            EC (IIPT,IRAGB) = EC (IIPT,IRAGB) + SCALE*V2RHOASIGMABB
            EC (IIPT,IRAGC) = EC (IIPT,IRAGC) + SCALE*V2RHOASIGMAAB
            EC (IIPT,IRBRB) = EC (IIPT,IRBRB) + SCALE*V2RHOB2
            EC (IIPT,IRBGA) = EC (IIPT,IRBGA) + SCALE*V2RHOBSIGMAAA
            EC (IIPT,IRBGB) = EC (IIPT,IRBGB) + SCALE*V2RHOBSIGMABB
            EC (IIPT,IRBGC) = EC (IIPT,IRBGC) + SCALE*V2RHOBSIGMAAB
            EC (IIPT,IGAGA) = EC (IIPT,IGAGA) + SCALE*V2SIGMAAA2
            EC (IIPT,IGAGB) = EC (IIPT,IGAGB) + SCALE*V2SIGMAAABB
            EC (IIPT,IGAGC) = EC (IIPT,IGAGC) + SCALE*V2SIGMAAAAB
            EC (IIPT,IGBGB) = EC (IIPT,IGBGB) + SCALE*V2SIGMABB2
            EC (IIPT,IGBGC) = EC (IIPT,IGBGC) + SCALE*V2SIGMAABBB
            EC (IIPT,IGCGC) = EC (IIPT,IGCGC) + SCALE*V2SIGMAAB2
C
         END IF
      END DO
C
      ECOR = ECOR + ECF
C
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFPBECF
      SUBROUTINE NDFPBECF(ECOR,EC0,EC,DRHOI_,XYZWGT,SCALE,RHOMIN,
     *                    ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      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 EC(ILENG,*)
      DIMENSION EC0(*)
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION XYZWGT(*)
      DATA ZERO/0.0D+00/
      DATA TWO /2.0D+00/
C
      IDERIV = NDER
C     IDERIV = NUMBER OF DERIVATIVES NEEDED
C              TDDFT ENERGIES REQUIRE NDER=2
C              TDDFT GRADIENTS REQUIRE NDER=3
C
      NPT    = 1
C     NPT    = NUMBER POINTS TO BE EVALUATED IN THE DFT REPOSITORY
C              ROUTINE
C
      ECF    = ZERO
C     ECF    = INITIALIZE THE ENERGY OF THE CORRELATION FUNCTION TO ZERO
C
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
         RHOA=DRHOI_(IPT,1,1)
         RHOB=DRHOI_(IPT,1,IRU)
C             IN THE CLOSED SHELL CASE RHOB =RHOA SINCE IRU=1
C             IN THE OPENED SHELL CASE RHOB!=RHOA SINCE IRU=2
         RHOT=RHOA + RHOB
C
         IF(RHOT/TWO.GT.RHOMIN) THEN
           DXA = DRHOI_(IPT,2,1)
           DYA = DRHOI_(IPT,3,1)
           DZA = DRHOI_(IPT,4,1)
           DXB = DRHOI_(IPT,2,IRU)
           DYB = DRHOI_(IPT,3,IRU)
           DZB = DRHOI_(IPT,4,IRU)
           GRDAA=DXA*DXA + DYA*DYA + DZA*DZA
           GRDBB=DXB*DXB + DYB*DYB + DZB*DZB
           GRDAB=DXA*DXB + DYA*DYB + DZA*DZB
C
           RHOA1=RHOA
           RHOB1=RHOB
           SIGMAAA1=GRDAA
           SIGMABB1=GRDBB
           SIGMAAB1=GRDAB
C
C          INITIALIZE FIRST DERIVATIVES
C
           VRHOA                = 0.00D+00
           VRHOB                = 0.00D+00
           VSIGMAAA             = 0.00D+00
           VSIGMABB             = 0.00D+00
           VSIGMAAB             = 0.00D+00
C
C          INITIALIZE SECOND DERIVATIVES
C
           V2RHOA2              = 0.00D+00
           V2RHOAB              = 0.00D+00
           V2RHOASIGMAAA        = 0.00D+00
           V2RHOASIGMABB        = 0.00D+00
           V2RHOASIGMAAB        = 0.00D+00
           V2RHOB2              = 0.00D+00
           V2RHOBSIGMAAA        = 0.00D+00
           V2RHOBSIGMABB        = 0.00D+00
           V2RHOBSIGMAAB        = 0.00D+00
           V2SIGMAAA2           = 0.00D+00
           V2SIGMAAABB          = 0.00D+00
           V2SIGMAAAAB          = 0.00D+00
           V2SIGMABB2           = 0.00D+00
           V2SIGMAABBB          = 0.00D+00
           V2SIGMAAB2           = 0.00D+00
C
           CALL UKS_C_PBE
     &    (IDERIV,NPT,RHOA1,RHOB1,SIGMAAA1,SIGMABB1,SIGMAAB1,
     &     ZK,VRHOA,VRHOB,VSIGMAAA,VSIGMABB,VSIGMAAB,
     &     V2RHOA2,V2RHOB2,V2RHOAB,
     &     V2RHOASIGMAAA,V2RHOASIGMAAB,V2RHOASIGMABB,
     &     V2RHOBSIGMABB,V2RHOBSIGMAAB,V2RHOBSIGMAAA,
     &     V2SIGMAAA2,V2SIGMAAAAB,V2SIGMAAABB,
     &     V2SIGMAAB2,V2SIGMAABBB,V2SIGMABB2)
C
            ECF             = ECF             + SCALE*ZK*XYZWGT(IPT)
            EC0(IIPT)       = EC0(IIPT)       + SCALE*ZK
C           FIRST DERIVATIVES
            EC (IIPT,IRA)   = EC (IIPT,IRA)   + SCALE*VRHOA
            EC (IIPT,IRB)   = EC (IIPT,IRB)   + SCALE*VRHOB
            EC (IIPT,IGA)   = EC (IIPT,IGA)   + SCALE*VSIGMAAA
            EC (IIPT,IGB)   = EC (IIPT,IGB)   + SCALE*VSIGMABB
            EC (IIPT,IGC)   = EC (IIPT,IGC)   + SCALE*VSIGMAAB
C           SECOND DERIVATIVES
            EC (IIPT,IRARA) = EC (IIPT,IRARA) + SCALE*V2RHOA2
            EC (IIPT,IRARB) = EC (IIPT,IRARB) + SCALE*V2RHOAB
            EC (IIPT,IRAGA) = EC (IIPT,IRAGA) + SCALE*V2RHOASIGMAAA
            EC (IIPT,IRAGB) = EC (IIPT,IRAGB) + SCALE*V2RHOASIGMABB
            EC (IIPT,IRAGC) = EC (IIPT,IRAGC) + SCALE*V2RHOASIGMAAB
            EC (IIPT,IRBRB) = EC (IIPT,IRBRB) + SCALE*V2RHOB2
            EC (IIPT,IRBGA) = EC (IIPT,IRBGA) + SCALE*V2RHOBSIGMAAA
            EC (IIPT,IRBGB) = EC (IIPT,IRBGB) + SCALE*V2RHOBSIGMABB
            EC (IIPT,IRBGC) = EC (IIPT,IRBGC) + SCALE*V2RHOBSIGMAAB
            EC (IIPT,IGAGA) = EC (IIPT,IGAGA) + SCALE*V2SIGMAAA2
            EC (IIPT,IGAGB) = EC (IIPT,IGAGB) + SCALE*V2SIGMAAABB
            EC (IIPT,IGAGC) = EC (IIPT,IGAGC) + SCALE*V2SIGMAAAAB
            EC (IIPT,IGBGB) = EC (IIPT,IGBGB) + SCALE*V2SIGMABB2
            EC (IIPT,IGBGC) = EC (IIPT,IGBGC) + SCALE*V2SIGMAABBB
            EC (IIPT,IGCGC) = EC (IIPT,IGCGC) + SCALE*V2SIGMAAB2
C
         END IF
      END DO
C
      ECOR = ECOR + ECF
C
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFVWN1
      SUBROUTINE NDFVWN1(ECOR,EC0,EC,DRHOI_,XYZWGT,SCALE,RHOMIN,
     *                   ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      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 EC(ILENG,*)
      DIMENSION EC0(*)
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION XYZWGT(*)
      DATA ZERO/0.0D+00/
      DATA TWO /2.0D+00/
C
      IDERIV = NDER
C     IDERIV = NUMBER OF DERIVATIVES NEEDED
C              TDDFT ENERGIES REQUIRE NDER=2
C              TDDFT GRADIENTS REQUIRE NDER=3
C
      NPT    = 1
C     NPT    = NUMBER POINTS TO BE EVALUATED IN THE DFT REPOSITORY
C              ROUTINE
C
      ECF    = ZERO
C     ECF    = INITIALIZE THE ENERGY OF THE CORRELATION FUNCTION TO ZERO
C
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
         RHOA=DRHOI_(IPT,1,1)
         RHOB=DRHOI_(IPT,1,IRU)
C             IN THE CLOSED SHELL CASE RHOB =RHOA SINCE IRU=1
C             IN THE OPENED SHELL CASE RHOB!=RHOA SINCE IRU=2
         RHOT=RHOA + RHOB
C
         IF(RHOT/TWO.GT.RHOMIN) THEN
C
           RHOA1=RHOA
           RHOB1=RHOB
C
C          INITIALIZE FIRST DERIVATIVES
C
           VRHOA                = 0.00D+00
           VRHOB                = 0.00D+00
C
C          INITIALIZE SECOND DERIVATIVES
C
           V2RHOA2              = 0.00D+00
           V2RHOAB              = 0.00D+00
           V2RHOB2              = 0.00D+00
C
           CALL UKS_C_VWN1RPA(IDERIV,NPT,RHOA1,RHOB1,ZK,
     &                        VRHOA,VRHOB,V2RHOA2,V2RHOB2,V2RHOAB)
C
            ECF             = ECF             + SCALE*ZK*XYZWGT(IPT)
            EC0(IIPT)       = EC0(IIPT)       + SCALE*ZK
C           FIRST DERIVATIVES
            EC (IIPT,IRA)   = EC (IIPT,IRA)   + SCALE*VRHOA
            EC (IIPT,IRB)   = EC (IIPT,IRB)   + SCALE*VRHOB
C           SECOND DERIVATIVES
            EC (IIPT,IRARA) = EC (IIPT,IRARA) + SCALE*V2RHOA2
            EC (IIPT,IRARB) = EC (IIPT,IRARB) + SCALE*V2RHOAB
            EC (IIPT,IRBRB) = EC (IIPT,IRBRB) + SCALE*V2RHOB2
C
         END IF
      END DO
C
      ECOR = ECOR + ECF
C
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFPW91LDA
      SUBROUTINE NDFPW91LDA(ECOR,EC0,EC,DRHOI_,XYZWGT,SCALE,RHOMIN,
     *                      ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      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 EC(ILENG,*)
      DIMENSION EC0(*)
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION XYZWGT(*)
      DATA ZERO/0.0D+00/
      DATA TWO /2.0D+00/
C
      IDERIV = NDER
C     IDERIV = NUMBER OF DERIVATIVES NEEDED
C              TDDFT ENERGIES REQUIRE NDER=2
C              TDDFT GRADIENTS REQUIRE NDER=3
C
      NPT    = 1
C     NPT    = NUMBER POINTS TO BE EVALUATED IN THE DFT REPOSITORY
C              ROUTINE
C
      ECF    = ZERO
C     ECF    = INITIALIZE THE ENERGY OF THE CORRELATION FUNCTION TO ZERO
C
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
         RHOA=DRHOI_(IPT,1,1)
         RHOB=DRHOI_(IPT,1,IRU)
C             IN THE CLOSED SHELL CASE RHOB =RHOA SINCE IRU=1
C             IN THE OPENED SHELL CASE RHOB!=RHOA SINCE IRU=2
         RHOT=RHOA + RHOB
C
         IF(RHOT/TWO.GT.RHOMIN) THEN
C
           RHOA1=RHOA
           RHOB1=RHOB
C
C          INITIALIZE FIRST DERIVATIVES
C
           VRHOA                = 0.00D+00
           VRHOB                = 0.00D+00
C
C          INITIALIZE SECOND DERIVATIVES
C
           V2RHOA2              = 0.00D+00
           V2RHOAB              = 0.00D+00
           V2RHOB2              = 0.00D+00
C
           CALL UKS_C_PW92(IDERIV,NPT,RHOA1,RHOB1,ZK,
     &                     VRHOA,VRHOB,V2RHOA2,V2RHOB2,V2RHOAB)
C
            ECF             = ECF             + SCALE*ZK*XYZWGT(IPT)
            EC0(IIPT)       = EC0(IIPT)       + SCALE*ZK
C           FIRST DERIVATIVES
            EC (IIPT,IRA)   = EC (IIPT,IRA)   + SCALE*VRHOA
            EC (IIPT,IRB)   = EC (IIPT,IRB)   + SCALE*VRHOB
C           SECOND DERIVATIVES
            EC (IIPT,IRARA) = EC (IIPT,IRARA) + SCALE*V2RHOA2
            EC (IIPT,IRARB) = EC (IIPT,IRARB) + SCALE*V2RHOAB
            EC (IIPT,IRBRB) = EC (IIPT,IRBRB) + SCALE*V2RHOB2
C
         END IF
      END DO
C
      ECOR = ECOR + ECF
C
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFP86CF
      SUBROUTINE NDFP86CF(ECOR,EC0,EC,DRHOI_,XYZWGT,SCALE,RHOMIN,
     *                    ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      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 EC(ILENG,*)
      DIMENSION EC0(*)
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION XYZWGT(*)
      DATA ZERO/0.0D+00/
      DATA TWO /2.0D+00/
C
      IDERIV = NDER
C     IDERIV = NUMBER OF DERIVATIVES NEEDED
C              TDDFT ENERGIES REQUIRE NDER=2
C              TDDFT GRADIENTS REQUIRE NDER=3
C
      NPT    = 1
C     NPT    = NUMBER POINTS TO BE EVALUATED IN THE DFT REPOSITORY
C              ROUTINE
C
      ECF    = ZERO
C     ECF    = INITIALIZE THE ENERGY OF THE CORRELATION FUNCTION TO ZERO
C
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
         RHOA=DRHOI_(IPT,1,1)
         RHOB=DRHOI_(IPT,1,IRU)
C             IN THE CLOSED SHELL CASE RHOB =RHOA SINCE IRU=1
C             IN THE OPENED SHELL CASE RHOB!=RHOA SINCE IRU=2
         RHOT=RHOA + RHOB
C
         IF(RHOT/TWO.GT.RHOMIN) THEN
           DXA = DRHOI_(IPT,2,1)
           DYA = DRHOI_(IPT,3,1)
           DZA = DRHOI_(IPT,4,1)
           DXB = DRHOI_(IPT,2,IRU)
           DYB = DRHOI_(IPT,3,IRU)
           DZB = DRHOI_(IPT,4,IRU)
           GRDAA=DXA*DXA + DYA*DYA + DZA*DZA
           GRDBB=DXB*DXB + DYB*DYB + DZB*DZB
           GRDAB=DXA*DXB + DYA*DYB + DZA*DZB
C
           RHOA1=RHOA
           RHOB1=RHOB
           SIGMAAA1=GRDAA
           SIGMABB1=GRDBB
           SIGMAAB1=GRDAB
C
C          INITIALIZE FIRST DERIVATIVES
C
           VRHOA                = 0.00D+00
           VRHOB                = 0.00D+00
           VSIGMAAA             = 0.00D+00
           VSIGMABB             = 0.00D+00
           VSIGMAAB             = 0.00D+00
C
C          INITIALIZE SECOND DERIVATIVES
C
           V2RHOA2              = 0.00D+00
           V2RHOAB              = 0.00D+00
           V2RHOASIGMAAA        = 0.00D+00
           V2RHOASIGMABB        = 0.00D+00
           V2RHOASIGMAAB        = 0.00D+00
           V2RHOB2              = 0.00D+00
           V2RHOBSIGMAAA        = 0.00D+00
           V2RHOBSIGMABB        = 0.00D+00
           V2RHOBSIGMAAB        = 0.00D+00
           V2SIGMAAA2           = 0.00D+00
           V2SIGMAAABB          = 0.00D+00
           V2SIGMAAAAB          = 0.00D+00
           V2SIGMABB2           = 0.00D+00
           V2SIGMAABBB          = 0.00D+00
           V2SIGMAAB2           = 0.00D+00
C
           CALL UKS_C_P86
     &    (IDERIV,NPT,RHOA1,RHOB1,SIGMAAA1,SIGMABB1,
     &     SIGMAAB1,ZK,VRHOA,VRHOB,VSIGMAAA,VSIGMABB,VSIGMAAB,
     &     V2RHOA2,V2RHOB2,V2RHOAB,
     &     V2RHOASIGMAAA,V2RHOASIGMAAB,V2RHOASIGMABB,
     &     V2RHOBSIGMABB,V2RHOBSIGMAAB,V2RHOBSIGMAAA,
     &     V2SIGMAAA2,V2SIGMAAAAB,V2SIGMAAABB,
     &     V2SIGMAAB2,V2SIGMAABBB,V2SIGMABB2)
C
            ECF             = ECF             + SCALE*ZK*XYZWGT(IPT)
            EC0(IIPT)       = EC0(IIPT)       + SCALE*ZK
C           FIRST DERIVATIVES
            EC (IIPT,IRA)   = EC (IIPT,IRA)   + SCALE*VRHOA
            EC (IIPT,IRB)   = EC (IIPT,IRB)   + SCALE*VRHOB
            EC (IIPT,IGA)   = EC (IIPT,IGA)   + SCALE*VSIGMAAA
            EC (IIPT,IGB)   = EC (IIPT,IGB)   + SCALE*VSIGMABB
            EC (IIPT,IGC)   = EC (IIPT,IGC)   + SCALE*VSIGMAAB
C           SECOND DERIVATIVES
            EC (IIPT,IRARA) = EC (IIPT,IRARA) + SCALE*V2RHOA2
            EC (IIPT,IRARB) = EC (IIPT,IRARB) + SCALE*V2RHOAB
            EC (IIPT,IRAGA) = EC (IIPT,IRAGA) + SCALE*V2RHOASIGMAAA
            EC (IIPT,IRAGB) = EC (IIPT,IRAGB) + SCALE*V2RHOASIGMABB
            EC (IIPT,IRAGC) = EC (IIPT,IRAGC) + SCALE*V2RHOASIGMAAB
            EC (IIPT,IRBRB) = EC (IIPT,IRBRB) + SCALE*V2RHOB2
            EC (IIPT,IRBGA) = EC (IIPT,IRBGA) + SCALE*V2RHOBSIGMAAA
            EC (IIPT,IRBGB) = EC (IIPT,IRBGB) + SCALE*V2RHOBSIGMABB
            EC (IIPT,IRBGC) = EC (IIPT,IRBGC) + SCALE*V2RHOBSIGMAAB
            EC (IIPT,IGAGA) = EC (IIPT,IGAGA) + SCALE*V2SIGMAAA2
            EC (IIPT,IGAGB) = EC (IIPT,IGAGB) + SCALE*V2SIGMAAABB
            EC (IIPT,IGAGC) = EC (IIPT,IGAGC) + SCALE*V2SIGMAAAAB
            EC (IIPT,IGBGB) = EC (IIPT,IGBGB) + SCALE*V2SIGMABB2
            EC (IIPT,IGBGC) = EC (IIPT,IGBGC) + SCALE*V2SIGMAABBB
            EC (IIPT,IGCGC) = EC (IIPT,IGCGC) + SCALE*V2SIGMAAB2
C
         END IF
      END DO
C
      ECOR = ECOR + ECF
C
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFPZ81CF
      SUBROUTINE NDFPZ81CF(ECOR,EC0,EC,DRHOI_,XYZWGT,SCALE,RHOMIN,
     *                     ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      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 EC(ILENG,*)
      DIMENSION EC0(*)
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION XYZWGT(*)
      DATA ZERO/0.0D+00/
      DATA TWO /2.0D+00/
C
      IDERIV = NDER
C     IDERIV = NUMBER OF DERIVATIVES NEEDED
C              TDDFT ENERGIES REQUIRE NDER=2
C              TDDFT GRADIENTS REQUIRE NDER=3
C
      NPT    = 1
C     NPT    = NUMBER POINTS TO BE EVALUATED IN THE DFT REPOSITORY
C              ROUTINE
C
      ECF    = ZERO
C     ECF    = INITIALIZE THE ENERGY OF THE CORRELATION FUNCTION TO ZERO
C
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
         RHOA=DRHOI_(IPT,1,1)
         RHOB=DRHOI_(IPT,1,IRU)
C             IN THE CLOSED SHELL CASE RHOB =RHOA SINCE IRU=1
C             IN THE OPENED SHELL CASE RHOB!=RHOA SINCE IRU=2
         RHOT=RHOA + RHOB
C
         IF(RHOT/TWO.GT.RHOMIN) THEN
C
           RHOA1=RHOA
           RHOB1=RHOB
C
C          INITIALIZE FIRST DERIVATIVES
C
           VRHOA                = 0.00D+00
           VRHOB                = 0.00D+00
C
C          INITIALIZE SECOND DERIVATIVES
C
           V2RHOA2              = 0.00D+00
           V2RHOAB              = 0.00D+00
           V2RHOB2              = 0.00D+00
C
           CALL UKS_C_PZ81(IDERIV,NPT,RHOA1,RHOB1,ZK,VRHOA,VRHOB,
     *                     V2RHOA2,V2RHOB2,V2RHOAB)
C
            ECF             = ECF             + SCALE*ZK*XYZWGT(IPT)
            EC0(IIPT)       = EC0(IIPT)       + SCALE*ZK
C           FIRST DERIVATIVES
            EC (IIPT,IRA)   = EC (IIPT,IRA)   + SCALE*VRHOA
            EC (IIPT,IRB)   = EC (IIPT,IRB)   + SCALE*VRHOB
C           SECOND DERIVATIVES
            EC (IIPT,IRARA) = EC (IIPT,IRARA) + SCALE*V2RHOA2
            EC (IIPT,IRARB) = EC (IIPT,IRARB) + SCALE*V2RHOAB
            EC (IIPT,IRBRB) = EC (IIPT,IRBRB) + SCALE*V2RHOB2
C
         END IF
      END DO
C
      ECOR = ECOR + ECF
C
      RETURN
      END
C
C     NDFNWPW91CF INTERFACES THE UNDOCUMENTED PW91 CORRELATION
C     FUNCTIONAL UTILIZING THE NWCHEM CONSTANTS FOR PW91LDA
C*MODULE TDFFUN  *DECK NDFNWPW91CF
      SUBROUTINE NDFNWPW91CF(ECOR,EC0,EC,DRHOI_,XYZWGT,SCALE,RHOMIN,
     *                       ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      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 EC(ILENG,*)
      DIMENSION EC0(*)
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION XYZWGT(*)
      DATA ZERO/0.0D+00/
      DATA TWO /2.0D+00/
C
      IDERIV = NDER
C     IDERIV = NUMBER OF DERIVATIVES NEEDED
C              TDDFT ENERGIES REQUIRE NDER=2
C              TDDFT GRADIENTS REQUIRE NDER=3
C
      NPT    = 1
C     NPT    = NUMBER POINTS TO BE EVALUATED IN THE DFT REPOSITORY
C              ROUTINE
C
      ECF    = ZERO
C     ECF    = INITIALIZE THE ENERGY OF THE CORRELATION FUNCTION TO ZERO
C
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
         RHOA=DRHOI_(IPT,1,1)
         RHOB=DRHOI_(IPT,1,IRU)
C             IN THE CLOSED SHELL CASE RHOB =RHOA SINCE IRU=1
C             IN THE OPENED SHELL CASE RHOB!=RHOA SINCE IRU=2
         RHOT=RHOA + RHOB
C
         IF(RHOT/TWO.GT.RHOMIN) THEN
           DXA = DRHOI_(IPT,2,1)
           DYA = DRHOI_(IPT,3,1)
           DZA = DRHOI_(IPT,4,1)
           DXB = DRHOI_(IPT,2,IRU)
           DYB = DRHOI_(IPT,3,IRU)
           DZB = DRHOI_(IPT,4,IRU)
           GRDAA=DXA*DXA + DYA*DYA + DZA*DZA
           GRDBB=DXB*DXB + DYB*DYB + DZB*DZB
           GRDAB=DXA*DXB + DYA*DYB + DZA*DZB
C
           RHOA1=RHOA
           RHOB1=RHOB
           SIGMAAA1=GRDAA
           SIGMABB1=GRDBB
           SIGMAAB1=GRDAB
C
C          INITIALIZE FIRST DERIVATIVES
C
           VRHOA                = 0.00D+00
           VRHOB                = 0.00D+00
           VSIGMAAA             = 0.00D+00
           VSIGMABB             = 0.00D+00
           VSIGMAAB             = 0.00D+00
C
C          INITIALIZE SECOND DERIVATIVES
C
           V2RHOA2              = 0.00D+00
           V2RHOAB              = 0.00D+00
           V2RHOASIGMAAA        = 0.00D+00
           V2RHOASIGMABB        = 0.00D+00
           V2RHOASIGMAAB        = 0.00D+00
           V2RHOB2              = 0.00D+00
           V2RHOBSIGMAAA        = 0.00D+00
           V2RHOBSIGMABB        = 0.00D+00
           V2RHOBSIGMAAB        = 0.00D+00
           V2SIGMAAA2           = 0.00D+00
           V2SIGMAAABB          = 0.00D+00
           V2SIGMAAAAB          = 0.00D+00
           V2SIGMABB2           = 0.00D+00
           V2SIGMAABBB          = 0.00D+00
           V2SIGMAAB2           = 0.00D+00
C
           CALL UKS_PW91C_NWCHEM
     & (IDERIV,NPT,RHOA1,RHOB1,SIGMAAA1,SIGMABB1,SIGMAAB1,
     &  ZK,VRHOA,VRHOB,VSIGMAAA,VSIGMABB,VSIGMAAB,
     &  V2RHOA2,V2RHOB2,V2RHOAB,
     &  V2RHOASIGMAAA,V2RHOASIGMAAB,V2RHOASIGMABB,
     &  V2RHOBSIGMABB,V2RHOBSIGMAAB,V2RHOBSIGMAAA,
     &  V2SIGMAAA2,V2SIGMAAAAB,V2SIGMAAABB,
     &  V2SIGMAAB2,V2SIGMAABBB,V2SIGMABB2)
C
            ECF             = ECF             + SCALE*ZK*XYZWGT(IPT)
            EC0(IIPT)       = EC0(IIPT)       + SCALE*ZK
C           FIRST DERIVATIVES
            EC (IIPT,IRA)   = EC (IIPT,IRA)   + SCALE*VRHOA
            EC (IIPT,IRB)   = EC (IIPT,IRB)   + SCALE*VRHOB
            EC (IIPT,IGA)   = EC (IIPT,IGA)   + SCALE*VSIGMAAA
            EC (IIPT,IGB)   = EC (IIPT,IGB)   + SCALE*VSIGMABB
            EC (IIPT,IGC)   = EC (IIPT,IGC)   + SCALE*VSIGMAAB
C           SECOND DERIVATIVES
            EC (IIPT,IRARA) = EC (IIPT,IRARA) + SCALE*V2RHOA2
            EC (IIPT,IRARB) = EC (IIPT,IRARB) + SCALE*V2RHOAB
            EC (IIPT,IRAGA) = EC (IIPT,IRAGA) + SCALE*V2RHOASIGMAAA
            EC (IIPT,IRAGB) = EC (IIPT,IRAGB) + SCALE*V2RHOASIGMABB
            EC (IIPT,IRAGC) = EC (IIPT,IRAGC) + SCALE*V2RHOASIGMAAB
            EC (IIPT,IRBRB) = EC (IIPT,IRBRB) + SCALE*V2RHOB2
            EC (IIPT,IRBGA) = EC (IIPT,IRBGA) + SCALE*V2RHOBSIGMAAA
            EC (IIPT,IRBGB) = EC (IIPT,IRBGB) + SCALE*V2RHOBSIGMABB
            EC (IIPT,IRBGC) = EC (IIPT,IRBGC) + SCALE*V2RHOBSIGMAAB
            EC (IIPT,IGAGA) = EC (IIPT,IGAGA) + SCALE*V2SIGMAAA2
            EC (IIPT,IGAGB) = EC (IIPT,IGAGB) + SCALE*V2SIGMAAABB
            EC (IIPT,IGAGC) = EC (IIPT,IGAGC) + SCALE*V2SIGMAAAAB
            EC (IIPT,IGBGB) = EC (IIPT,IGBGB) + SCALE*V2SIGMABB2
            EC (IIPT,IGBGC) = EC (IIPT,IGBGC) + SCALE*V2SIGMAABBB
            EC (IIPT,IGCGC) = EC (IIPT,IGCGC) + SCALE*V2SIGMAAB2
C
         END IF
      END DO
C
      ECOR = ECOR + ECF
C
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFPW91XF
      SUBROUTINE NDFPW91XF(EXC,EX0,EX,DRHOI_,XYZWGT,SCALE,RHOMIN,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      COMMON /XDERIX/ NXDIM(3),KRA,KRB,KGA,KGB,KRARA,KRAGA,KGAGA,
     *                KRBRB,KRBGB,KGBGB,KRARARA,KRARAGA,KRAGAGA,KGAGAGA,
     *                KRBRBRB,KRBRBGB,KRBGBGB,KGBGBGB
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION EX0(*)
      DIMENSION EX(ILENG,*)
      DIMENSION XYZWGT(*)
      DATA ZERO/0.0D+00/
      DATA TWO/2.0D+00/
C
C
      IDERIV = NDER
C     IDERIV = NUMBER OF DERIVATIVES NEEDED
C              TDDFT ENERGIES REQUIRE NDER=2
C              TDDFT GRADIENTS REQUIRE NDER=3
C
      NPT    = 1
C     NPT    = NUMBER POINTS TO BE EVALUATED IN THE DFT REPOSITORY
C              ROUTINE
C
      EXF    = ZERO
C     EXF    = INITIALIZE THE ENERGY OF THE CORRELATION FUNCTION TO ZERO
C
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
         RHOA=DRHOI_(IPT,1,1)
         RHOB=DRHOI_(IPT,1,IRU)
C             IN THE CLOSED SHELL CASE RHOB =RHOA SINCE IRU=1
C             IN THE OPENED SHELL CASE RHOB!=RHOA SINCE IRU=2
         RHOT=RHOA + RHOB
C
         IF(RHOT/TWO.GT.RHOMIN) THEN
           DXA = DRHOI_(IPT,2,1)
           DYA = DRHOI_(IPT,3,1)
           DZA = DRHOI_(IPT,4,1)
           DXB = DRHOI_(IPT,2,IRU)
           DYB = DRHOI_(IPT,3,IRU)
           DZB = DRHOI_(IPT,4,IRU)
           GRDAA=DXA*DXA + DYA*DYA + DZA*DZA
           GRDBB=DXB*DXB + DYB*DYB + DZB*DZB
C
           RHOA1=RHOA
           RHOB1=RHOB
           SIGMAAA1=GRDAA
           SIGMABB1=GRDBB
C
C          INITIALIZE FIRST DERIVATIVES
C
           VRHOA                = 0.00D+00
           VRHOB                = 0.00D+00
           VSIGMAAA             = 0.00D+00
           VSIGMABB             = 0.00D+00
C
C          INITIALIZE SECOND DERIVATIVES
C
           V2RHOA2              = 0.00D+00
           V2RHOAB              = 0.00D+00
           V2RHOASIGMAAA        = 0.00D+00
           V2RHOASIGMABB        = 0.00D+00
           V2RHOB2              = 0.00D+00
           V2RHOBSIGMAAA        = 0.00D+00
           V2RHOBSIGMABB        = 0.00D+00
           V2SIGMAAA2           = 0.00D+00
           V2SIGMABB2           = 0.00D+00
           V2SIGMAAB2           = 0.00D+00
C
           CALL UKS_X_PW91
     &    (IDERIV,NPT,RHOA1,RHOB1,SIGMAAA1,SIGMABB1,
     &     ZK,VRHOA,VRHOB,VSIGMAAA,VSIGMABB,VSIGMAAB,
     &     V2RHOA2,V2RHOB2,V2RHOAB,
     &     V2RHOASIGMAAA,V2RHOASIGMAAB,V2RHOASIGMABB,
     &     V2RHOBSIGMABB,V2RHOBSIGMAAB,V2RHOBSIGMAAA,
     &     V2SIGMAAA2,V2SIGMAAAAB,V2SIGMAAABB,
     &     V2SIGMAAB2,V2SIGMAABBB,V2SIGMABB2)
C
            EXF             = EXF             + SCALE*ZK*XYZWGT(IPT)
            EX0(IIPT)       = EX0(IIPT)       + SCALE*ZK
C           FIRST DERIVATIVES
            EX (IIPT,KRA)   = EX (IIPT,KRA)   + SCALE*VRHOA
            EX (IIPT,KRB)   = EX (IIPT,KRB)   + SCALE*VRHOB
            EX (IIPT,KGA)   = EX (IIPT,KGA)   + SCALE*VSIGMAAA
            EX (IIPT,KGB)   = EX (IIPT,KGB)   + SCALE*VSIGMABB
C           SECOND DERIVATIVES
            EX (IIPT,KRARA) = EX (IIPT,KRARA) + SCALE*V2RHOA2
            EX (IIPT,KRAGA) = EX (IIPT,KRAGA) + SCALE*V2RHOASIGMAAA
            EX (IIPT,KRBRB) = EX (IIPT,KRBRB) + SCALE*V2RHOB2
            EX (IIPT,KRBGB) = EX (IIPT,KRBGB) + SCALE*V2RHOBSIGMABB
            EX (IIPT,KGAGA) = EX (IIPT,KGAGA) + SCALE*V2SIGMAAA2
            EX (IIPT,KGBGB) = EX (IIPT,KGBGB) + SCALE*V2SIGMABB2
         END IF
      END DO
C
      EXC = EXC + EXF
C
      RETURN
      END
C*MODULE TDFFUN  *DECK NDFPBEXF
      SUBROUTINE NDFPBEXF(EXC,EX0,EX,DRHOI_,XYZWGT,SCALE,RHOMIN,
     *               ILENG,NPTGRD,NDER,IRU,IST,IEND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      COMMON /XDERIX/ NXDIM(3),KRA,KRB,KGA,KGB,KRARA,KRAGA,KGAGA,
     *                KRBRB,KRBGB,KGBGB,KRARARA,KRARAGA,KRAGAGA,KGAGAGA,
     *                KRBRBRB,KRBRBGB,KRBGBGB,KGBGBGB
      DIMENSION DRHOI_(NPTGRD,4,IRU)
      DIMENSION EX0(*)
      DIMENSION EX(ILENG,*)
      DIMENSION XYZWGT(*)
      DATA ZERO/0.0D+00/
      DATA TWO/2.0D+00/
C
      IDERIV = NDER
C     IDERIV = NUMBER OF DERIVATIVES NEEDED
C              TDDFT ENERGIES REQUIRE NDER=2
C              TDDFT GRADIENTS REQUIRE NDER=3
C
      NPT    = 1
C     NPT    = NUMBER POINTS TO BE EVALUATED IN THE DFT REPOSITORY
C              ROUTINE
C
      EXF    = ZERO
C     EXF    = INITIALIZE THE ENERGY OF THE CORRELATION FUNCTION TO ZERO
C
      DO IPT=IST,IEND
         IIPT=IPT-IST+1
         RHOA=DRHOI_(IPT,1,1)
         RHOB=DRHOI_(IPT,1,IRU)
C             IN THE CLOSED SHELL CASE RHOB =RHOA SINCE IRU=1
C             IN THE OPENED SHELL CASE RHOB!=RHOA SINCE IRU=2
         RHOT=RHOA + RHOB
C
         IF(RHOT/TWO.GT.RHOMIN) THEN
           DXA = DRHOI_(IPT,2,1)
           DYA = DRHOI_(IPT,3,1)
           DZA = DRHOI_(IPT,4,1)
           DXB = DRHOI_(IPT,2,IRU)
           DYB = DRHOI_(IPT,3,IRU)
           DZB = DRHOI_(IPT,4,IRU)
           GRDAA=DXA*DXA + DYA*DYA + DZA*DZA
           GRDBB=DXB*DXB + DYB*DYB + DZB*DZB
C
           RHOA1=RHOA
           RHOB1=RHOB
           SIGMAAA1=GRDAA
           SIGMABB1=GRDBB
C
C          INITIALIZE FIRST DERIVATIVES
C
           VRHOA                = 0.00D+00
           VRHOB                = 0.00D+00
           VSIGMAAA             = 0.00D+00
           VSIGMABB             = 0.00D+00
C
C          INITIALIZE SECOND DERIVATIVES
C
           V2RHOA2              = 0.00D+00
           V2RHOAB              = 0.00D+00
           V2RHOASIGMAAA        = 0.00D+00
           V2RHOASIGMABB        = 0.00D+00
           V2RHOB2              = 0.00D+00
           V2RHOBSIGMAAA        = 0.00D+00
           V2RHOBSIGMABB        = 0.00D+00
           V2SIGMAAA2           = 0.00D+00
           V2SIGMABB2           = 0.00D+00
           V2SIGMAAB2           = 0.00D+00
C
           CALL UKS_X_PBE
     &    (IDERIV,NPT,RHOA1,RHOB1,SIGMAAA1,SIGMABB1,
     &     ZK,VRHOA,VRHOB,VSIGMAAA,VSIGMABB,VSIGMAAB,
     &     V2RHOA2,V2RHOB2,V2RHOAB,
     &     V2RHOASIGMAAA,V2RHOASIGMAAB,V2RHOASIGMABB,
     &     V2RHOBSIGMABB,V2RHOBSIGMAAB,V2RHOBSIGMAAA,
     &     V2SIGMAAA2,V2SIGMAAAAB,V2SIGMAAABB,
     &     V2SIGMAAB2,V2SIGMAABBB,V2SIGMABB2)
C
            EXF             = EXF             + SCALE*ZK*XYZWGT(IPT)
            EX0(IIPT)       = EX0(IIPT)       + SCALE*ZK
C           FIRST DERIVATIVES
            EX (IIPT,KRA)   = EX (IIPT,KRA)   + SCALE*VRHOA
            EX (IIPT,KRB)   = EX (IIPT,KRB)   + SCALE*VRHOB
            EX (IIPT,KGA)   = EX (IIPT,KGA)   + SCALE*VSIGMAAA
            EX (IIPT,KGB)   = EX (IIPT,KGB)   + SCALE*VSIGMABB
C           SECOND DERIVATIVES
            EX (IIPT,KRARA) = EX (IIPT,KRARA) + SCALE*V2RHOA2
            EX (IIPT,KRAGA) = EX (IIPT,KRAGA) + SCALE*V2RHOASIGMAAA
            EX (IIPT,KRBRB) = EX (IIPT,KRBRB) + SCALE*V2RHOB2
            EX (IIPT,KRBGB) = EX (IIPT,KRBGB) + SCALE*V2RHOBSIGMABB
            EX (IIPT,KGAGA) = EX (IIPT,KGAGA) + SCALE*V2SIGMAAA2
            EX (IIPT,KGBGB) = EX (IIPT,KGBGB) + SCALE*V2SIGMABB2
         END IF
      END DO
C
      EXC = EXC + EXF
C
      RETURN
      END
C
