C  4 MAR 08 - SS  - NEW FILE FOR META-GGA FUNCTIONALS
C=======================================================================
C     THE FOLLOWING FUNCTIONAL CODES AND THEIR VARIANTS: M05, M06, VS98
C     WERE ORIGINALLY WRITTEN BY TRUHLAR THEORETICAL AND COMPUTATIONAL
C     CHEMISTRY GROUP AT THE UNIVERSITY OF MINNESOTA.
C
C     THEIR CODE WAS INTEGRATED INTO GAMESS WITH THEIR PERMISSION BY
C     SAROM SOK AT IOWA STATE UNIVERSITY.
C=======================================================================
C*MODULE DFTXCC *DECK M05CF
      SUBROUTINE M05CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                 GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                 FTOTWT,ECF,
     >                 VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA,
     >                 VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB,
     >                 TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IJZY)
C     [A] ZHAO, Y.; SCHULTZ, N. E.; TRUHLAR, D. G. J. CHEM. PHYS. 2005,
C         123, 161103;
C     NOTE THAT IN THIS COMMUNICATION WE INTERCHANGED CCAB,I AND CCSS,I
C     IN TABLE 1.
C     [B] ZHAO, Y.; SCHULTZ, N. E.; TRUHLAR, D. G. J. CHEM. THEORY
C         COMPUT. 2006, IN PRESS.
      IMPLICIT NONE
C
C     CLOSED-SHELL OPEN-SHELL LOGICALS
C
      LOGICAL RHFGVB,UROHF
C
C     EXCHANGE ENERGY
C
      DOUBLE PRECISION ECF
C
C     CHARGE DENSITY
C
      DOUBLE PRECISION ROA,ROB
C
C     CHARGE DENSITY GRADIENT AND COMPONENTS
C
      DOUBLE PRECISION GRDAA,GRDBB
      DOUBLE PRECISION GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB
C
C     KINETIC ENERGY DENSITY COMPONENTS
C
      DOUBLE PRECISION TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB
C
C     QUADRATURE WEIGHTS
C
      DOUBLE PRECISION FTOTWT
C
C     THE XC POTENTIAL
C
      DOUBLE PRECISION VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA
      DOUBLE PRECISION VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB
C
C     PARAMETER CONTROL
C
      INTEGER IJZY
C
C     WORKING VARIABLES
C
      DOUBLE PRECISION F1, F2, F3, F4,COPP
      DATA COPP/0.0031D0/
      DATA F1/1.0D0/,F2/2.0D0/,F3/3.0D0/,F4/4.0D0/
C
C    CALL TO THE M05CSS SUBROUTINE
C
      DOUBLE PRECISION PA,GAA,TA,FA,FPA,FGA,FTA,EUA,CHIA,EUPA
      DOUBLE PRECISION CHIAP,CHIAG
      DOUBLE PRECISION PB,GBB,TB,FB,FPB,FGB,FTB,EUB,CHIB,EUPB
      DOUBLE PRECISION CHIBP,CHIBG
C
      DOUBLE PRECISION SOP, SOPP0, SOPP1,SOPP2, SOPP3, SOPP4
      DOUBLE PRECISION PI, PI34, F13
      DOUBLE PRECISION RS,RSP,ZETA,DZDA,DZDB,POTLC,DLDS,DLDZ,D2LDSS
      DOUBLE PRECISION D2LDSZ,D2LDZZ
      DOUBLE PRECISION P, EUEG, U, W
      DOUBLE PRECISION DUDCHIA,DUDCHIB,DUDPA,DUDPB,DUDGA,DUDGB
      DOUBLE PRECISION DWDU,DWDPA,DWDPB, DWDGA,DWDGB,EUEGPA,EUEGPB
C
C     SPIN-RESTRICTED
C
      SOP=1.0D0
      IF (IJZY.EQ.1) THEN
C     PARAMETERS FOR M05 CORRELATION
         SOPP0 =  1.00000D0
         SOPP1 =  3.78569D0
         SOPP2 =-14.15261D0
         SOPP3 = -7.46589D0
         SOPP4 = 17.94491D0
      ELSEIF (IJZY.EQ.2) THEN
C     PARAMETERS FOR M05-2X CORRELATION
         SOPP0 =  1.00000D0
         SOPP1 =  1.09297D0
         SOPP2 = -3.79171D0
         SOPP3 =  2.82810D0
         SOPP4 =-10.58909D0
      ENDIF
C
      PI   = F4*ATAN(F1)
      PI34 = F3 / (F4*PI)
      F13  = F1 / F3
C
      IF (RHFGVB) THEN
C
C     GET THE DENSITY, GRADIENT, AND TAU FOR THE ALPHA SPIN FROM THE
C     TOTAL
C
        PA  = ROA
        GAA = GRDAA
C
C     IN THE M05CSS SUBROUTINE, WE USE 2*TA AS THE TAU, SO WE DO NOT
C     DIVIDE THE TAU BY 2 HERE
C
        TA = TAUXA + TAUYA + TAUZA
C
        CALL M05CSS(PA,GAA,TA,FA,FPA,FGA,FTA,EUA,
     >              CHIA,EUPA,CHIAP,CHIAG,IJZY)
        PB    = PA
        GBB   = GAA
        TB    = TA
        FB    = FA
        FPB   = FPA
        FGB   = FGA
        FTB   = FTA
        EUB   = EUA
        CHIB  = CHIA
        EUPB  = EUPA
        CHIBP = CHIAP
        CHIBG = CHIAG
C
        ECF = ECF + 2.0D+00*FA*FTOTWT
C
        VXCA1 = VXCA1 + FPA
C
        DUMAX = DUMAX + 2.0D+00*(FGA*GRADXA)
        DUMAY = DUMAY + 2.0D+00*(FGA*GRADYA)
        DUMAZ = DUMAZ + 2.0D+00*(FGA*GRADZA)
C
        DMXCA = DMXCA + FTA
      END IF
C
C     SPIN-UNRESTRICTED
C
      IF (UROHF) THEN
C
C     ALPHA
C
        PA  = ROA
        GAA = GRDAA
C
C     IN THE M05CSS SUBROUTINE, WE USE 2*TA AS THE TAU
C
        TA = TAUXA + TAUYA + TAUZA
C
        CALL M05CSS(PA,GAA,TA,FA,FPA,FGA,FTA,EUA,
     >              CHIA,EUPA,CHIAP,CHIAG,IJZY)
C
        ECF = ECF + FA*FTOTWT
C
        VXCA1 = VXCA1 + FPA
C
        DUMAX = DUMAX + 2.0D0*(FGA*GRADXA)
        DUMAY = DUMAY + 2.0D0*(FGA*GRADYA)
        DUMAZ = DUMAZ + 2.0D0*(FGA*GRADZA)
C
        DMXCA = DMXCA + FTA
C
C     BETA
C
        PB = ROB
        GBB = GRDBB
C
C     IN THE M05CSS SUBROUTINE, WE USE 2*TA AS THE TAU
C
        TB = TAUXB + TAUYB + TAUZB
C
        CALL M05CSS(PB,GBB,TB,FB,FPB,FGB,FTB,EUB,
     >              CHIB,EUPB,CHIBP,CHIBG,IJZY)
C
        ECF = ECF + FB*FTOTWT
C
        VXCB1 = VXCB1 + FPB
C
        DUMBX = DUMBX + 2.0D0*(FGB*GRADXB)
        DUMBY = DUMBY + 2.0D0*(FGB*GRADYB)
        DUMBZ = DUMBZ + 2.0D0*(FGB*GRADZB)
C
        DMXCB = DMXCB + FTB
      ENDIF
C
      P = PA + PB
C
      RS   =(PI34/P) ** F13
      RSP  =-RS/(F3*P)
      ZETA =(PA-PB)/P
      DZDA =(F1-ZETA)/P
      DZDB =(-F1-ZETA)/P
C
      CALL LSDAC(RS,ZETA,POTLC,DLDS,DLDZ,D2LDSS,D2LDSZ,
     >           D2LDZZ)
C
      EUEG = P*POTLC - EUA - EUB
      U    = COPP*(CHIA+CHIB)/(F1 + COPP*(CHIA+CHIB))
      W    = SOPP0+U*(SOPP1+U*(SOPP2+U*(SOPP3+U*SOPP4)))
C
      ECF = ECF + SOP*EUEG*W*FTOTWT
C
      DUDCHIA =COPP/(F1 + COPP*(CHIA+CHIB))**2
      DUDCHIB =COPP/(F1 + COPP*(CHIA+CHIB))**2
C
      DUDPA = DUDCHIA*CHIAP
      DUDPB = DUDCHIB*CHIBP
      DUDGA = DUDCHIA*CHIAG
      DUDGB = DUDCHIB*CHIBG
      DWDU  = SOPP1+U*(F2*SOPP2+U*(F3*SOPP3+U*F4*SOPP4))
      DWDPA = DWDU*DUDPA
      DWDPB = DWDU*DUDPB
      DWDGA = DWDU*DUDGA
      DWDGB = DWDU*DUDGB
C
      EUEGPA = POTLC + P*DLDS*RSP + P*DLDZ*DZDA - EUPA
      EUEGPB = POTLC + P*DLDS*RSP + P*DLDZ*DZDB - EUPB
C
        IF (RHFGVB) THEN
          VXCA1 = VXCA1 + SOP*(EUEGPA*W + EUEG*DWDPA)
C
          DUMAX = DUMAX + 2.0D+00*SOP*(EUEG*DWDGA)*GRADXA
          DUMAY = DUMAY + 2.0D+00*SOP*(EUEG*DWDGA)*GRADYA
          DUMAZ = DUMAZ + 2.0D+00*SOP*(EUEG*DWDGA)*GRADZA
        ELSE
          VXCA1 = VXCA1 + SOP*(EUEGPA*W + EUEG*DWDPA)
          VXCB1 = VXCB1 + SOP*(EUEGPB*W + EUEG*DWDPB)
C
          DUMAX = DUMAX + 2.0D+00*SOP*EUEG*DWDGA*GRADXA
          DUMAY = DUMAY + 2.0D+00*SOP*EUEG*DWDGA*GRADYA
          DUMAZ = DUMAZ + 2.0D+00*SOP*EUEG*DWDGA*GRADZA
C
          DUMBX = DUMBX + 2.0D+00*SOP*EUEG*DWDGB*GRADXB
          DUMBY = DUMBY + 2.0D+00*SOP*EUEG*DWDGB*GRADYB
          DUMBZ = DUMBZ + 2.0D+00*SOP*EUEG*DWDGB*GRADZB
        ENDIF
      END
C
C*MODULE DFTXCC *DECK M05CSS
      SUBROUTINE M05CSS(PX,GX,TX,F,FP,FG,FT,EUEG,CHI,EUEGP,
     >                  CHIP,CHIG,IJZY)
      IMPLICIT NONE
C
C     COMPUTE THE SAME-SPIN PART OF THE M05 CORRELATION FUNCTIONAL FOR
C     ONE GRID POINT AND ONE SPIN-CASE.
C
      INTEGER IJZY
      DOUBLE PRECISION PX, GX, TX, F, FP, FG, FT
      DOUBLE PRECISION EUEG, CHI, EUEGP, CHIP, CHIG
      DOUBLE PRECISION PT25, F1, F2, F3, F4, F8
      DOUBLE PRECISION SS, SSS0,SSS1, SSS2, SSS3, SSS4, CSS
      DOUBLE PRECISION PI, PI34, F13, F83
      DOUBLE PRECISION RS, D, FSCC, RSP, DFSCCP, DFSCCG
      DOUBLE PRECISION E, W, U, DFSCCT, DUDCHI, DWDU, DWDP, DWDG
      DOUBLE PRECISION D2LDSS,D2LDSZ,D2LDZZ,POTLC,DLDS,DLDZ
C
      DATA PT25/0.25D0/, F1/1.0D0/, F2/2.0D0/, F3/3.0D0/,
     >     F4/4.0D0/, F8/8.0D0/, CSS/0.06D0/
C
      SS = 1.0D+00
C
      IF (IJZY.EQ.1) THEN
C     PARAMETERS FOR M05 CORRELATION
         SSS0 =   1.00000D0
         SSS1 =   3.77344D0
         SSS2 = -26.04463D0
         SSS3 =  30.69913D0
         SSS4 =  -9.22695D0
      ELSEIF (IJZY.EQ.2) THEN
C     PARAMETERS FOR M05-2X CORRELATION
         SSS0 =   1.00000D0
         SSS1 =  -3.05430D0
         SSS2 =   7.61854D0
         SSS3 =   1.47665D0
         SSS4 = -11.92365D0
      ENDIF
C
      PI   = F4*ATAN(F1)
      PI34 = F3 / (F4*PI)
      F13  = F1 / F3
      F83  = F8 / F3
C
C
      RS = (PI34/PX) ** F13
C
      CALL LSDAC(RS,F1,POTLC,DLDS,DLDZ,D2LDSS,D2LDSZ,D2LDZZ)
C
      EUEG = PX*POTLC
C
      D = TX - PT25*GX/PX
C
      CHI = GX/PX**F83
C
      U = CSS*CHI/(F1 + CSS*CHI)
      W = SSS0+U*(SSS1+U*(SSS2+U*(SSS3+U*SSS4)))
C
      FSCC = D/TX
C
      E = FSCC*W*EUEG
      F = E*SS
C
      RSP = -RS/(F3*PX)
C
      CHIG = F1/PX**F83
      CHIP = -F83*CHI/PX
C
      DFSCCP = PT25*GX/(TX*PX**2)
      DFSCCG =-PT25/(TX*PX)
      DFSCCT = PT25*GX/(PX*TX**2)
      DUDCHI = CSS/((F1+CSS*CHI)**2)
C
      DWDU = SSS1+U*(F2*SSS2+U*(F3*SSS3+U*F4*SSS4))
      DWDP = DWDU*DUDCHI*CHIP
      DWDG = DWDU*DUDCHI*CHIG
C
      EUEGP = POTLC + PX*DLDS*RSP
C
      FP = SS*(DFSCCP*W*EUEG
     >   + FSCC*DWDP*EUEG
     >   + FSCC*W*EUEGP)
C
      FG = SS*(DFSCCG*W*EUEG
     >   + FSCC*DWDG*EUEG)
C
      FT = SS*(DFSCCT*W*EUEG)
      RETURN
      END
C
C     PERDEW 91 LOCAL CORRELATION FUNCTIONAL AT ONE GRID POINT
C
C*MODULE DFTXCC *DECK LSDCA
      SUBROUTINE LSDAC(RS,ZETA,POTLC,DLDS,DLDZ,D2LDSS,D2LDSZ,D2LDZZ)
      IMPLICIT NONE
C
C     EVALUATE THE PERDEW 91 LOCAL CORRELATION FUNCTIONAL AND ITS
C     DERIVATIVES AT ONE POINT.
C
      DOUBLE PRECISION RS
      DOUBLE PRECISION POTLC,DLDS, DLDZ
      DOUBLE PRECISION EPS0C(6), EPS1C(6), EPSC(6)
      DOUBLE PRECISION F1, F2, F3, F4, F6, F8, F9, F12
      DOUBLE PRECISION GAMMAI,ZETA,FZETA,DFZDZ,D2FZDZ
      DOUBLE PRECISION EU,DEUDRS,D2UDRS
      DOUBLE PRECISION EP,DEPDRS,D2PDRS
      DOUBLE PRECISION ALPHAM,DAMDRS,D2ADRS
      DOUBLE PRECISION GZ, HZ, DGZ, DHZ, D2GZ, D2HZ
C
      DOUBLE PRECISION PI,THIRD,Z2,Z3,Z4,D2LDSS,D2LDSZ,D2LDZZ,FZZI
C
      DATA EPS0C/0.03109070D0,0.21370D0, 7.5957D0,3.5876D0,1.6382D0,
     >           0.49294D0/
      DATA EPS1C/0.01554535D0,0.20548D0,14.1189D0,6.1977D0,3.3662D0,
     >           0.62517D0/
      DATA EPSC /0.01688690D0,0.11125D0,10.3570D0,3.6231D0,0.88026D0,
     >           0.49671D0/
      DATA F1/1.0D0/, F2/2.0D0/, F3/3.0D0/, F4/4.0D0/, F6/6.0D0/,
     >     F8/8.0D0/, F9/9.0D0/, F12/12.0D0/
C
      PI    = F4*ATAN(F1)
      THIRD = F1 / F3
C
      EPSC(1) = F1 / (F6*PI*PI)
      FZZI    = F9*(F2**THIRD-F1) / F4
      GAMMAI  = F1 / (F2*F2**THIRD-F2)
C
      CALL EVFZET(GAMMAI,ZETA,FZETA,DFZDZ,D2FZDZ)
C
      CALL EVPWLC(EPS0C(1),EPS0C(2),EPS0C(3),EPS0C(4),EPS0C(5),
     >            EPS0C(6), RS,EU,DEUDRS,D2UDRS)
C
      CALL EVPWLC(EPS1C(1),EPS1C(2),EPS1C(3),EPS1C(4),EPS1C(5),
     >            EPS1C(6),RS,EP,DEPDRS,D2PDRS)
C
      CALL EVPWLC(EPSC(1),EPSC(2),EPSC(3),EPSC(4),EPSC(5),
     >            EPSC(6),RS,ALPHAM,DAMDRS,D2ADRS)
C
      Z2 = ZETA*ZETA
      Z3 = ZETA*Z2
      Z4 = ZETA*Z3
      GZ = FZETA*Z4
      HZ = FZZI*(FZETA-GZ)
C
      POTLC = EU*(F1-GZ) + EP*GZ - ALPHAM*HZ
C
      DLDS = DEUDRS*(F1-GZ) + DEPDRS*GZ - DAMDRS*HZ
      DGZ  = DFZDZ*Z4 + F4*FZETA*Z3
      DHZ  = FZZI*(DFZDZ-DGZ)
      DLDZ = (EP-EU)*DGZ - ALPHAM*DHZ
      D2GZ = D2FZDZ*Z4 + F8*Z3*DFZDZ + F12*FZETA*Z2
      D2HZ = FZZI*(D2FZDZ-D2GZ)
C
      D2LDSS = D2UDRS*(F1-GZ) + D2PDRS*GZ - D2ADRS*HZ
      D2LDSZ = (DEPDRS-DEUDRS)*DGZ - DAMDRS*DHZ
      D2LDZZ = (EP-EU)*D2GZ - ALPHAM*D2HZ
      RETURN
      END
C
C     F(ZETA)
C
C*MODULE DFTXCC *DECK EVFZET
      SUBROUTINE EVFZET(S,ZETA,FZETA,DFZDZ,D2FZDZ)
      IMPLICIT REAL*8(A-H,O-Z)
C
C     EVALUATE F(ZETA) AND ITS DERIVATIVES FOR LSDAC.
C
      DOUBLE PRECISION SMALL
      DOUBLE PRECISION S, ZETA, FZETA,DFZDZ,D2FZDZ
      DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, NINE
      DOUBLE PRECISION OMZ, OPZ, OMZ3, OPZ3
      DOUBLE PRECISION F13, F43, F49
      DATA ZERO/0.0D0/, ONE/1.0D0/, TWO/2.0D0/, THREE/3.0D0/,
     >     FOUR/4.0D0/, NINE/9.0D0/
C
      SMALL  = 1.0D-14
      FZETA  =-TWO
      DFZDZ  = ZERO
      D2FZDZ = ZERO
C
      OMZ  = ONE - ZETA
      OPZ  = ONE + ZETA
      F13  = ONE / THREE
      F43  = FOUR / THREE
      F49  = FOUR / NINE
C
      IF(OMZ.GT.SMALL) THEN
        OMZ3   = OMZ ** F13
        FZETA  = FZETA + OMZ*OMZ3
        DFZDZ  = DFZDZ - OMZ3
        D2FZDZ = D2FZDZ + OMZ3/OMZ
      ENDIF
C
      IF(OPZ.GT.SMALL) THEN
        OPZ3   = OPZ ** F13
        FZETA  = FZETA + OPZ*OPZ3
        DFZDZ  = DFZDZ + OPZ3
        D2FZDZ = D2FZDZ + OPZ3/OPZ
      ENDIF
C
      FZETA  = FZETA * S
      DFZDZ  = DFZDZ * F43 * S
      D2FZDZ = D2FZDZ * F49 * S
      RETURN
      END
C
C     PW91 LOCAL CORRELATION
C
C*MODULE DFTXCC *DECK EVPWLC
      SUBROUTINE EVPWLC(A,A1,B1,B2,B3,B4,RS,V,DVDRS,D2VDRS)
      IMPLICIT NONE
C
C     EVALUATE THE INTERPOLATION FUNCTION FOR PW91 LOCAL CORRELATION.
C
      DOUBLE PRECISION A,A1,B1,B2,B3,B4,RS,V,DVDRS,D2VDRS
      DOUBLE PRECISION F1,F2, F3, F4
      DOUBLE PRECISION Q0, RS12, RS32,Q1,Q2
      DOUBLE PRECISION DQ0DRS,DQ1DRS,DQ2DRS
      DOUBLE PRECISION D2Q1DS, D2Q2DS
      DATA F1/1.0D0/, F2/2.0D0/, F3/3.0D0/, F4/4.0D0/
C
      Q0   =-F2*A*(F1+A1*RS)
      RS12 = SQRT(RS)
      RS32 = RS*RS12
      Q1   = F2*A*(B1*RS12+B2*RS+B3*RS32+B4*RS*RS)
      Q2   = LOG(F1+F1/Q1)
      V    = Q0*Q2
C
      DQ0DRS =-F2*A*A1
      DQ1DRS = A*(B1/RS12+F2*B2+F3*B3*RS12+F4*B4*RS)
      DQ2DRS =-DQ1DRS/(Q1+Q1**2)
      DVDRS  = DQ0DRS*Q2 + Q0*DQ2DRS
C
      D2Q1DS = A*(F3*B3/(RS12*F2)-B1/(RS32*F2)+F4*B4)
      D2Q2DS = (F2*Q1+F1)*(DQ1DRS/(Q1+Q1**2))**2 - D2Q1DS/(Q1+Q1**2)
      D2VDRS = F2*DQ0DRS*DQ2DRS + Q0*D2Q2DS
      RETURN
      END
C*MODULE DFTXCC *DECK M05XF
      SUBROUTINE M05XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                 GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                 FTOTWT,XGRD,
     >                 VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA,
     >                 VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB,
     >                 TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IJZY)
C     [A] ZHAO, Y.; SCHULTZ, N. E.; TRUHLAR, D. G. J. CHEM. PHYS. 2005,
C         123, 161103;
C     NOTE THAT IN THIS COMMUNICATION WE INTERCHANGED CCAB,I AND CCSS,I
C     IN TABLE 1.
C     [B] ZHAO, Y.; SCHULTZ, N. E.; TRUHLAR, D. G. J. CHEM. THEORY
C         COMPUT. 2006, IN PRESS.
C
      IMPLICIT NONE
C
      INTEGER IM05
      COMMON /CM05  / CM05XF, IM05
C
C     CLOSED-SHELL OPEN-SHELL LOGICALS
C
      LOGICAL RHFGVB,UROHF
C
C     EXCHANGE ENERGY
C
      DOUBLE PRECISION XGRD
C
C     CHARGE DENSITY
C
      DOUBLE PRECISION ROA,ROB
C
C     CHARGE DENSITY GRADIENT AND COMPONENTS
C
      DOUBLE PRECISION GRDAA,GRDBB
      DOUBLE PRECISION GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB
C
C     KINETIC ENERGY DENSITY COMPONENTS
C
      DOUBLE PRECISION TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB
C
C     QUADRATURE WEIGHTS
C
      DOUBLE PRECISION FTOTWT
C
C     THE XC POTENTIAL
C
      DOUBLE PRECISION VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA
      DOUBLE PRECISION VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB
C
C     PARAMETER CONTROL
C
      INTEGER IJZY
C
C     SCALE FACTORS FOR LOCAL AND NON-LOCAL CONTRIBUTIONS.
C
      DOUBLE PRECISION CM05XF
C
C     WORKING VARIABLES
C
      DOUBLE PRECISION PI
C
      DOUBLE PRECISION AT1, AT2, AT3, AT4, AT5, AT6, AT7, AT8, AT9
      DOUBLE PRECISION AT, AT10, AT11, C1, C2, FL, FNL
      DOUBLE PRECISION RRHO, RHO43, RHO13, RHOO, RHO53
      DOUBLE PRECISION GAMMA2, GAMMA
      DOUBLE PRECISION TAUUEG, TSIG, WSIG, W1, W2, W3, W4, W5, W6
      DOUBLE PRECISION W7, W8, W9, W10, W11, FSIG
C
      DOUBLE PRECISION TAUN,TAUU
      DOUBLE PRECISION F23, F53, F13
      DOUBLE PRECISION F4O3
      DOUBLE PRECISION ONE, TWO, THREE, FOUR, FIVE, SIX, SEVEN, EIGHT
      DOUBLE PRECISION NINE, F10, F11
      DOUBLE PRECISION CS, AX, P32, S, S2, EN, ED, E, DE, DEN, DED
C
C      FUNCTIONAL DERIVATIVES
C
      DOUBLE PRECISION DFDW, DWDT, DTDR, DTDTAU, DGGADR, DFDR
      DOUBLE PRECISION DFDTAU, DGGADG
C
      PARAMETER (PI = 3.1415926535897932384626433832795D0)
C
      PARAMETER (F13=1.D0/3.D0)
      PARAMETER (F4O3=4.D0/3.D0)
      PARAMETER (F23=2.0D0/3.D0, F53=5.D0/3.D0)
      PARAMETER (ONE=1.0D0, TWO=2.0D0, THREE=3.0D0, FOUR=4.0D0,
     >           FIVE=5.0D0,SIX=6.0D0, SEVEN=7.0D0,
     >           EIGHT=8.0D0, NINE=9.0D0,F10=10.D0, F11=11.D0)
C
      FL  = CM05XF
      FNL = CM05XF
C
      IF (IJZY.EQ.1) THEN
C     PARAMETERS FOR M05
        AT1 =  0.08151D0
        AT2 = -0.43956D0
        AT3 = -3.22422D0
        AT4 =  2.01819D0
        AT5 =  8.79431D0
        AT6 = -0.00295D0
        AT7 =  9.82029D0
        AT8 = -4.82351D0
        AT9 =-48.17574D0
        AT10=  3.64802D0
        AT11= 34.02248D0
      ELSEIF (IJZY.EQ.2) THEN
C     PARAMETERS FOR M05-2X
        AT1 = -0.56833D0
        AT2 = -1.30057D0
        AT3 =  5.50070D0
        AT4 =  9.06402D0
        AT5 =-32.21075D0
        AT6 =-23.73298D0
        AT7 = 70.22996D0
        AT8 = 29.88614D0
        AT9 =-60.25778D0
        AT10=-13.22205D0
        AT11= 15.23694D0
      ENDIF
C
      AT = 1.0D0
      C1 = 0.2195149727645171D0
      C2 = C1/0.804D0
C
      CS  = 0.5D0/(3.0D0*PI*PI)**F13
      P32 = (3.D0*PI**2)**F23
C
      AX = (-0.75D0)*(3.0D0/PI)**F13
C
      IF (RHFGVB) THEN
        RHOO  = ROA+ROB
        RHO43 = RHOO**F4O3
        RRHO  = 1D0/RHOO
        RHO13 = RHO43*RRHO
        RHO53 = RHOO**F53
C
        TAUN = TAUXA + TAUYA + TAUZA
C
        TAUU   = TAUN
        TAUUEG = 0.3D0*P32*RHO53
        TSIG   = TAUUEG/TAUU
        WSIG   =(TSIG-ONE)/(TSIG+ONE)
        W1     = WSIG
        W2     = WSIG*W1
        W3     = WSIG*W2
        W4     = WSIG*W3
        W5     = WSIG*W4
        W6     = WSIG*W5
        W7     = WSIG*W6
        W8     = WSIG*W7
        W9     = WSIG*W8
        W10    = WSIG*W9
        W11    = WSIG*W10
        FSIG   = AT*(AT1*W1+ AT2*W2 + AT3*W3
     >         + AT4*W4 + AT5*W5 + AT6*W6 + AT7*W7
     >         + AT8*W8 + AT9*W9 + AT10*W10 + AT11*W11)
C
        GAMMA2 = 4.0D+00*GRDAA
        GAMMA  = SQRT(GAMMA2)
        S      = CS*GAMMA/RHO43
        S2     = S*S
        EN     = C1*S2
        ED     = ONE + C2*S2
        E      = EN/ED
C
        XGRD = XGRD + RHO43*AX*(FL+FNL*E)*(ONE+FSIG)*FTOTWT
C
C     FUNCTIONAL DERIVATIVES
C
        DEN    = TWO*C1*S
        DED    = TWO*C2*S
        DE     =(DEN*ED-EN*DED)/(ED*ED)
        DFDW   = AT*(AT1 + TWO*AT2*W1 + THREE*AT3*W2
     >         + FOUR*AT4*W3 + FIVE*AT5*W4
     >         + SIX*AT6*W5 + SEVEN*AT7*W6
     >         + EIGHT*AT8*W7 + NINE*AT9*W8
     >         + F10*AT10*W9 + F11*AT11*W10)
        DWDT   = TWO/((ONE + TSIG)**2)
        DTDR   =(0.5D0*P32*RHO13*RHO13)/TAUU
        DTDTAU =-TAUUEG/TAUU**2
        DGGADR = F4O3*RHO13*AX*(FL+FNL*(E-S*DE))
        DFDR   = DFDW*DWDT*DTDR
        DFDTAU = DFDW*DWDT*DTDTAU
        DGGADG =(FNL*DE*S/(TWO*GAMMA2))
C
        VXCA1 = VXCA1 + DGGADR*(ONE+FSIG)+(FL+FNL*E)*AX*RHO43*DFDR
C
        DUMAX = DUMAX + 2.0D+00*(TWO*DGGADG*AX*RHO43*(ONE+FSIG)*GRADXA)
        DUMAY = DUMAY + 2.0D+00*(TWO*DGGADG*AX*RHO43*(ONE+FSIG)*GRADYA)
        DUMAZ = DUMAZ + 2.0D+00*(TWO*DGGADG*AX*RHO43*(ONE+FSIG)*GRADZA)
C
        DMXCA = DMXCA + 0.5D0*RHO43*AX*(FL+FNL*E)*DFDTAU
      END IF
C
      IF (UROHF) THEN
C
C     USE SPIN DENSITY FUNCTIONAL THEORY IE N-->2N
C     EX=(1/2)EX[2*ALPHA] + (1/2)EX[2*BETA]
C
C     ALPHA
C
        RHOO  = TWO*ROA
        RHO43 = RHOO**F4O3
        RRHO  = 1.0D0/RHOO
        RHO13 = RHO43*RRHO
        RHO53 = RHOO**F53
C
        TAUN   = (TAUXA + TAUYA + TAUZA)*(ONE/TWO)
        TAUU   = TWO*TAUN
        TAUUEG = 0.3D0*P32*RHO53
        TSIG   = TAUUEG/TAUU
        WSIG   =(TSIG-ONE)/(TSIG+ONE)
        W1     = WSIG
        W2     = WSIG*W1
        W3     = WSIG*W2
        W4     = WSIG*W3
        W5     = WSIG*W4
        W6     = WSIG*W5
        W7     = WSIG*W6
        W8     = WSIG*W7
        W9     = WSIG*W8
        W10    = WSIG*W9
        W11    = WSIG*W10
        FSIG   = AT*(AT1*W1+ AT2*W2 + AT3*W3
     >         + AT4*W4 + AT5*W5 + AT6*W6 + AT7*W7
     >         + AT8*W8 + AT9*W9 + AT10*W10 + AT11*W11)
C
        GAMMA2 = GRDAA
        GAMMA2 = FOUR*GAMMA2
        GAMMA  = SQRT(GAMMA2)
C
        S      = CS*GAMMA/RHO43
        S2     = S*S
        EN     = C1*S2
        ED     = ONE + C2*S2
        E      = EN/ED
C
        XGRD = XGRD + RHO43*AX*(FL+FNL*E)*(ONE+FSIG)*0.5D0*FTOTWT
C
C     FUNCTIONAL DERIVATIVES
C
        DEN    = TWO*C1*S
        DED    = TWO*C2*S
        DE     =(DEN*ED-EN*DED)/(ED*ED)
        DFDW   = AT*(AT1 + TWO*AT2*W1 + THREE*AT3*W2
     >         + FOUR*AT4*W3 + FIVE*AT5*W4
     >         + SIX*AT6*W5 + SEVEN*AT7*W6
     >         + EIGHT*AT8*W7 + NINE*AT9*W8
     >         + F10*AT10*W9 + F11*AT11*W10)
        DWDT   = TWO/((ONE + TSIG)**2)
        DTDR   = TWO*(0.5D0*P32*RHO13*RHO13)/TAUU
        DTDTAU =-TWO*TAUUEG/TAUU**2
        DGGADR = TWO*F4O3*RHO13*AX*(FL+FNL*(E-S*DE))
        DFDR   = DFDW*DWDT*DTDR
        DFDTAU = DFDW*DWDT*DTDTAU
        DGGADG = FOUR*(FNL*DE*S/(TWO*GAMMA2))
C
        VXCA1 = VXCA1 + (DGGADR*(ONE+FSIG)+(FL+FNL*E)*AX*RHO43*DFDR)*
     >                   0.5D0
C
        DUMAX = DUMAX + 2.0D0*
     >                 (DGGADG*AX*RHO43*(ONE+FSIG)*0.5D0*GRADXA)
        DUMAY = DUMAY + 2.0D0*
     >                 (DGGADG*AX*RHO43*(ONE+FSIG)*0.5D0*GRADYA)
        DUMAZ = DUMAZ + 2.0D0*
     >                 (DGGADG*AX*RHO43*(ONE+FSIG)*0.5D0*GRADZA)
C
        DMXCA = DMXCA + 0.5D0*0.5D0*RHO43*AX*(FL+FNL*E)*DFDTAU
C
C     BETA
C
        RHOO  = TWO*ROB
        RHO43 = RHOO**F4O3
        RRHO  = 1.0D0/RHOO
        RHO13 = RHO43*RRHO
        RHO53 = RHOO**F53
C
        TAUN   = (TAUXB + TAUYB + TAUZB)*(ONE/TWO)
        TAUU   = TWO*TAUN
        TAUUEG = 0.3D0*P32*RHO53
        TSIG   = TAUUEG/TAUU
        WSIG   =(TSIG-ONE)/(TSIG+ONE)
        W1     = WSIG
        W2     = WSIG*W1
        W3     = WSIG*W2
        W4     = WSIG*W3
        W5     = WSIG*W4
        W6     = WSIG*W5
        W7     = WSIG*W6
        W8     = WSIG*W7
        W9     = WSIG*W8
        W10    = WSIG*W9
        W11    = WSIG*W10
        FSIG   = AT*(AT1*W1+ AT2*W2 + AT3*W3
     >         + AT4*W4 + AT5*W5 + AT6*W6 + AT7*W7
     >         + AT8*W8 + AT9*W9 + AT10*W10 + AT11*W11)
C
        GAMMA2 = GRDBB
        GAMMA2 = FOUR*GAMMA2
        GAMMA  = SQRT(GAMMA2)
        S      = CS*GAMMA/RHO43
        S2     = S*S
        EN     = C1*S2
        ED     = ONE + C2*S2
        E      = EN/ED
C
        XGRD = XGRD + RHO43*AX*(FL+FNL*E)*(ONE+FSIG)*0.5D0*FTOTWT
C
C     FUNCTIONAL DERIVATIVES
C
        DEN    = TWO*C1*S
        DED    = TWO*C2*S
        DE     = (DEN*ED-EN*DED)/(ED*ED)
        DFDW   = AT*(AT1 + TWO*AT2*W1 + THREE*AT3*W2
     >         + FOUR*AT4*W3 + FIVE*AT5*W4
     >         + SIX*AT6*W5 + SEVEN*AT7*W6
     >         + EIGHT*AT8*W7 + NINE*AT9*W8
     >         + F10*AT10*W9 + F11*AT11*W10)
        DWDT   = TWO/((ONE + TSIG)**2)
        DTDR   = TWO*(0.5D0*P32*RHO13*RHO13)/TAUU
        DTDTAU = -TWO*TAUUEG/TAUU**2
        DGGADR = TWO*F4O3*RHO13*AX*(FL+FNL*(E-S*DE))
        DFDR   = DFDW*DWDT*DTDR
        DFDTAU = DFDW*DWDT*DTDTAU
        DGGADG = FOUR*(FNL*DE*S/(TWO*GAMMA2))
C
        VXCB1 = VXCB1 + (DGGADR*(ONE+FSIG)+(FL+FNL*E)*AX*RHO43*DFDR)*
     >                 0.5D0
C
        DUMBX = DUMBX + 2.0D0*
     >                 (DGGADG*AX*RHO43*(ONE+FSIG)*0.5D0*GRADXB)
        DUMBY = DUMBY + 2.0D0*
     >                 (DGGADG*AX*RHO43*(ONE+FSIG)*0.5D0*GRADYB)
        DUMBZ = DUMBZ + 2.0D0*
     >                 (DGGADG*AX*RHO43*(ONE+FSIG)*0.5D0*GRADZB)
C
        DMXCB = DMXCB + 0.5D0*0.5D0*RHO43*AX*(FL+FNL*E)*DFDTAU
C
      ENDIF
      RETURN
      END
C*MODULE DFTXCC *DECK VS98CF
      SUBROUTINE VS98CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                  GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                  FTOTWT,ECF,
     >                  VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA,
     >                  VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB,
     >                  TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IJZY)
C
C     REFERENCE
C     [A] T. V. VOORHIS & G. E. SCUSERIA, J. CHEM. PHYS. 109, 400 (1998)
C     [B] Y. ZHAO & D. G. TRUHLAR, J. CHEM. PHYS. 125, 194101 (2006)
C
      IMPLICIT NONE
C
C     THRESHOLD PARAMETERS
C
      DOUBLE PRECISION F1, F3, F4, GAB
C
      DATA F1/1.0D+00/,F3/3.0D+00/,F4/4.0D+00/,GAB/0.00304966D+00/
C
C     CORRELATION ENERGY
C
      DOUBLE PRECISION ECF
C
C     CHARGE DENSITY
C
      DOUBLE PRECISION ROA,ROB
C
C     CHARGE DENSITY GRADIENT
C
      DOUBLE PRECISION GRDAA,GRDBB
      DOUBLE PRECISION GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB
C
C     KINETIC ENERGY DENSITY
C
      DOUBLE PRECISION TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB
C
C     QUADRATURE WEIGHTS
C
      DOUBLE PRECISION FTOTWT
C
C     THE XC POTENTIAL
C
      DOUBLE PRECISION VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA
      DOUBLE PRECISION VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB
C
      INTEGER IJZY
C
C     CALL TO THE VS98CSS SUBROUTINE
C
      DOUBLE PRECISION PA,GAA,TA,FA,FPA,FGA,FTA,EUA,CHIA,EUPA,
     >                 CHIAP,CHIAG,ZA,ZAP,ZAT
      DOUBLE PRECISION PB,GBB,TB,FB,FPB,FGB,FTB,EUB,CHIB,EUPB,
     >                 CHIBP,CHIBG,ZB,ZBP,ZBT
C
      DOUBLE PRECISION PI, F13, PI34, POTLC,RS,RSP,
     >                 ZETA,DZDA,DZDB,DLDS,DLDZ,D2LDSS,D2LDSZ,D2LDZZ
C
      DOUBLE PRECISION P, EUEG, ZAB, XAB, KAB, XK, ZK
      DOUBLE PRECISION DGDX,DGDZ,DGDPA,DGDGA,DGDTA,DGDPB,DGDGB,DGDTB
      DOUBLE PRECISION EUEGPA,EUEGPB,GCAB
      DOUBLE PRECISION R7, R8, R9, R10, R11, R12
C
      LOGICAL RHFGVB,UROHF
C
C     PARAMETERS FOR VS98
      IF (IJZY.EQ.1) THEN
              R7  = 7.035010D-01
              R8  = 7.694574D-03
              R9  = 5.152765D-02
              R10 = 3.394308D-05
              R11 =-1.269420D-03
              R12 = 1.296118D-03
C     PARAMETERS FOR M06-L
      ELSEIF (IJZY.EQ.2) THEN
              R7  = 3.957626D-01
              R8  =-5.614546D-01
              R9  = 1.403963D-02
              R10 = 9.831442D-04
              R11 =-3.577176D-03
              R12 = 0.000000D+00
C     PARAMETERS FOR M06-HF
      ELSEIF (IJZY.EQ.3) THEN
              R7  =-6.746338D-01
              R8  =-1.534002D-01
              R9  =-9.021521D-02
              R10 =-1.292037D-03
              R11 =-2.352983D-04
              R12 = 0.000000D+00
C     PARAMETERS FOR M06
      ELSEIF (IJZY.EQ.4) THEN
              R7  =-2.741539D+00
              R8  =-6.720113D-01
              R9  =-7.932688D-02
              R10 = 1.918681D-03
              R11 =-2.032902D-03
              R12 = 0.000000D+00
C     PARAMETERS FOR M06-2X
      ELSEIF (IJZY.EQ.5) THEN
              R7  = 1.166404D-01
              R8  =-9.120847D-02
              R9  =-6.726189D-02
              R10 = 6.720580D-05
              R11 = 8.448011D-04
              R12 = 0.000000D+00
      ENDIF
C
      PI   = F4*ATAN(F1)
      PI34 = F3 / (F4*PI)
      F13  = F1 / F3
C
      IF (RHFGVB) THEN
C
C     GET DENSITY, GRADIENT, AND TAU FOR THE ALPHA SPIN FROM THE TOTAL
C
        PA  = ROA
        GAA = GRDAA
C
C     IN THE BC95CSS SUBROUTINE, WE USE 2*TA AS THE TAU, SO WE DO NOT
C     DIVIDE THE TAU BY 2 HERE
C
        TA = TAUXA + TAUYA + TAUZA
C
        CALL VS98SS(PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
     >              CHIA,EUPA,CHIAP,CHIAG,ZAP,ZAT,IJZY)
C
        PB    = PA
        GBB   = GAA
        TB    = TA
        FB    = FA
        FPB   = FPA
        FGB   = FGA
        FTB   = FTA
        EUB   = EUA
        ZB    = ZA
        CHIB  = CHIA
        EUPB  = EUPA
        CHIBP = CHIAP
        CHIBG = CHIAG
        ZBP   = ZAP
        ZBT   = ZAT
C
        ECF = ECF + 2.0D+00*FA*FTOTWT
C
        VXCA1 = VXCA1 + FPA
C
        DUMAX = DUMAX + 2.0D+00*(FGA*GRADXA)
        DUMAY = DUMAY + 2.0D+00*(FGA*GRADYA)
        DUMAZ = DUMAZ + 2.0D+00*(FGA*GRADZA)
C
        DMXCA = DMXCA + FTA
      END IF
C
      IF (UROHF) THEN
C
C     ALPHA
C
        PA  = ROA
        GAA = GRDAA
C
C     IN THE BC95CSS SUBROUTINE, WE USE 2*TA AS THE TAU
C
        TA = TAUXA + TAUYA + TAUZA
C
        CALL VS98SS(PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
     >              CHIA,EUPA,CHIAP,CHIAG,ZAP,ZAT,IJZY)
C
        ECF = ECF + FA*FTOTWT
C
        VXCA1 = VXCA1 + FPA
C
        DUMAX = DUMAX + 2.0D+00*(FGA*GRADXA)
        DUMAY = DUMAY + 2.0D+00*(FGA*GRADYA)
        DUMAZ = DUMAZ + 2.0D+00*(FGA*GRADZA)
C
        DMXCA = DMXCA + FTA
C
C     BETA
C
        PB  = ROB
        GBB = GRDBB
C
C     IN THE VS98SS SUBROUTINE, WE USE 2*TA AS THE TAU,
C
        TB = TAUXB + TAUYB + TAUZB
C
        CALL VS98SS(PB,GBB,TB,FB,FPB,FGB,FTB,EUB,ZB,
     >              CHIB,EUPB,CHIBP,CHIBG,ZBP,ZBT,IJZY)
C
        ECF = ECF + FB*FTOTWT
C
        VXCB1 = VXCB1 + FPB
C
        DUMBX = DUMBX + 2.0D+00*(FGB*GRADXB)
        DUMBY = DUMBY + 2.0D+00*(FGB*GRADYB)
        DUMBZ = DUMBZ + 2.0D+00*(FGB*GRADZB)
C
        DMXCB = DMXCB + FTB
      ENDIF
C
      P = ROA + ROB
C
      RS   =(PI34/P) ** F13
      RSP  =-RS/(F3*P)
      ZETA =(PA-PB)/P
      DZDA =(F1-ZETA)/P
      DZDB =(-F1-ZETA)/P
C
      CALL LSDAC(RS,ZETA,POTLC,DLDS,DLDZ,D2LDSS,D2LDSZ,
     >           D2LDZZ)
C
      EUEG = P*POTLC - EUA - EUB
      ZAB  = ZA + ZB
      XAB  = CHIA+CHIB
      KAB  = F1 + GAB*(XAB+ZAB)
      XK   = XAB/KAB
      ZK   = ZAB/KAB
C
      CALL GVT4(GCAB,DGDX,DGDZ,XK,ZK,KAB,GAB,GAB,R7,R8,R9,R10,R11,R12)
C
      ECF = ECF + GCAB*EUEG*FTOTWT
C
      DGDPA = DGDX*CHIAP + DGDZ*ZAP
      DGDGA = DGDX*CHIAG
      DGDTA = DGDZ*ZAT
      DGDPB = DGDX*CHIBP + DGDZ*ZBP
      DGDGB = DGDX*CHIBG
      DGDTB = DGDZ*ZBT
C
      EUEGPA = POTLC + P*DLDS*RSP + P*DLDZ*DZDA - EUPA
      EUEGPB = POTLC + P*DLDS*RSP + P*DLDZ*DZDB - EUPB
C
      IF (RHFGVB) THEN
        VXCA1 = VXCA1 + (EUEGPA*GCAB + EUEG*DGDPA)
C
        DUMAX = DUMAX + 2.0D+00*(EUEG*DGDGA*GRADXA)
        DUMAY = DUMAY + 2.0D+00*(EUEG*DGDGA*GRADYA)
        DUMAZ = DUMAZ + 2.0D+00*(EUEG*DGDGA*GRADZA)
C
        DMXCA = DMXCA + EUEG*DGDTA
      END IF
C
      IF (UROHF) THEN
        VXCA1 = VXCA1 + (EUEGPA*GCAB + EUEG*DGDPA)
        VXCB1 = VXCB1 + (EUEGPB*GCAB + EUEG*DGDPB)
C
        DUMAX = DUMAX + 2.0D+00*(EUEG*DGDGA*GRADXA)
        DUMAY = DUMAY + 2.0D+00*(EUEG*DGDGA*GRADYA)
        DUMAZ = DUMAZ + 2.0D+00*(EUEG*DGDGA*GRADZA)
C
        DUMBX = DUMBX + 2.0D+00*(EUEG*DGDGB*GRADXB)
        DUMBY = DUMBY + 2.0D+00*(EUEG*DGDGB*GRADYB)
        DUMBZ = DUMBZ + 2.0D+00*(EUEG*DGDGB*GRADZB)
C
        DMXCA = DMXCA + EUEG*DGDTA
        DMXCB = DMXCB + EUEG*DGDTB
      ENDIF
      END
C
C*MODULE DFTXCC *DECK VS98SS
      SUBROUTINE VS98SS(PX,GX,TX,F,FP,FG,FT,EUEG,Z,CHI,EUEGP,
     >                  CHIP,CHIG,ZP,ZT,IJZY)
      IMPLICIT NONE
C
C     COMPUTE THE SAME-SPIN PART OF THE VS98 CORRELATION FUNCTIONAL
C     FOR ONE GRID POINT AND ONE SPIN-CASE.
C
      INTEGER IJZY
      DOUBLE PRECISION R13, R14, R15, R16, R17, R18
      DOUBLE PRECISION PX, GX, TX, F, FP, FG, FT, Z, ZP, ZT
      DOUBLE PRECISION EUEG, CHI, EUEGP, CHIP, CHIG, CF, GCC
      DOUBLE PRECISION F1, F3, F4, F5, F8
      DOUBLE PRECISION PI, PI34, F13, F53, F83
      DOUBLE PRECISION RS, D, RSP, POTLC, DX, DZ, DGDP, DGDG, DGDT
      DOUBLE PRECISION E,DP, DG, DT, RHOO, RHO53, RHO83
      DOUBLE PRECISION KC, XK, ZK, GC, DGDX, DGDZ
      DOUBLE PRECISION D2LDSS, D2LDSZ, D2LDZZ, DLDS, DLDZ
C
      DATA F1/1.0D+00/, F3/3.0D+00/, F4/4.0D+00/, F5/5.0D+00/,
     >     F8/8.0D+00/, GCC/0.00515088D+00/, CF/9.115599720D+00/
C
C     PARAMETERS FOR VS98
      IF (IJZY.EQ.1) THEN
              R13=   3.270912D-01
              R14=  -3.228915D-02
              R15=  -2.942406D-02
              R16=   2.134222D-03
              R17=  -5.451559D-03
              R18=   1.577575D-02
C     PARAMETERS FOR M06-L
      ELSEIF (IJZY.EQ.2) THEN
              R13=   4.650534D-01
              R14=   1.617589D-01
              R15=   1.833657D-01
              R16=   4.692100D-04
              R17=  -4.990573D-03
              R18=   0.000000D+00
C     PARAMETERS FOR M06-HF
      ELSEIF (IJZY.EQ.3) THEN
              R13=   8.976746D-01
              R14=  -2.345830D-01
              R15=   2.368173D-01
              R16=  -9.913890D-04
              R17=  -1.146165D-02
              R18=   0.000000D+00
C     PARAMETERS FOR M06
      ELSEIF (IJZY.EQ.4) THEN
               R13=  4.905945D-01
               R14= -1.437348D-01
               R15=  2.357824D-01
               R16=  1.871015D-03
               R17= -3.788963D-03
               R18=  0.000000D+00
C     PARAMETERS FOR M06-2X
      ELSEIF (IJZY.EQ.5) THEN
              R13=  6.902145D-01
              R14=  9.847204D-02
              R15=  2.214797D-01
              R16= -1.968264D-03
              R17= -6.775479D-03
              R18=  0.000000D+00
      ENDIF
C
      PI    = F4*ATAN(F1)
      PI34  = F3 / (F4*PI)
      F13   = F1 / F3
      F53   = F5 / F3
      F83   = F8 / F3
      RHOO  = PX
      RHO53 = RHOO**F53
      RHO83 = RHO53*RHOO
C
      RS = (PI34/PX) ** F13
C
      CALL LSDAC(RS,F1,POTLC,DLDS,DLDZ,D2LDSS,D2LDSZ,D2LDZZ)
C
      EUEG = PX*POTLC
      CHI  = GX/RHO83
      Z    =(TX/RHO53) - CF
      KC   = F1 + GCC*(CHI + Z)
      XK   = CHI/KC
      ZK   = Z/KC
      D    = F1 - CHI/(F4*(Z + CF))
C
      CALL GVT4(GC,DGDX,DGDZ,XK,ZK,KC,GCC,GCC,R13,R14,R15,R16,R17,R18)
C
      E = D*EUEG*GC
      F = E
C
      RSP   =-RS/(F3*PX)
      CHIG  = F1/PX**F83
      CHIP  =-F83*CHI/PX
      ZP    =-F53 * TX/RHO83
      ZT    = F1/RHO53
      DZ    = CHI/(F4*(Z + CF)*(Z + CF))
      DX    =-F1/(F4*(Z + CF))
      DP    = DZ*ZP + DX*CHIP
      DG    = DX*CHIG
      DT    = DZ*ZT
      DGDP  = DGDX*CHIP + DGDZ*ZP
      DGDG  = DGDX*CHIG
      DGDT  = DGDZ*ZT
      EUEGP = POTLC + PX*DLDS*RSP
      FP    = DP*EUEG*GC + D*EUEGP*GC + D*EUEG*DGDP
      FG    = DG*EUEG*GC + D*EUEG*DGDG
      FT    = DT*EUEG*GC + D*EUEG*DGDT
      RETURN
      END
C
C*MODULE DFTXCC *DECK GVT4
      SUBROUTINE GVT4(G,DGDX,DGDZ,XK,ZK,K,C,CT,R1,R2,R3,R4,R5,R6)
      IMPLICIT NONE
C
C     EVALUATE THE GVT4 FORM IN THE VS98 FUNCTIONAL
C
C     SOME WORKING VARIABLES
C
      DOUBLE PRECISION G,DGDX,DGDZ,XK,ZK,K,C,CT,R1,R2,R3,R4,R5,R6
      DOUBLE PRECISION SXK,SZK,SK,SC,SCT,SR1,SR2,SR3,SR4,SR5,SR6,SK2
      DOUBLE PRECISION ONE, TWO, THREE
      DATA ONE/1.0D+00/, TWO/2.0D+00/, THREE/3.0D+00/
C
      SXK  = XK
      SZK  = ZK
      SK   = K
      SC   = C
      SCT  = CT
      SR1  = R1
      SR2  = R2
      SR3  = R3
      SR4  = R4
      SR5  = R5
      SR6  = R6
      SK2  = SK*SK
C
      G    =(SR1 + SR2*SXK + SR3*SZK
     >     + SR4*SXK*SXK
     >     + SR5*SZK*SXK
     >     + SR6*SZK*SZK)/SK
      DGDX =(-SR1*SC
     >     +  SR2*(ONE-TWO*SC*SXK)
     >     -  TWO*SR3*SZK*SC
     >     +  SR4*(TWO*SXK-THREE*SXK*SXK*SC)
     >     +  SR5*(SZK -THREE*SZK*SXK*SC)
     >     -  THREE*SR6*SZK*SZK*SC)/SK2
      DGDZ =(-SR1*SCT
     >     -  TWO*SR2*SXK*SCT
     >     +  SR3*(ONE-TWO*SZK*SCT)
     >     -  THREE*SR4*SXK*SXK*SCT
     >     +  SR5*(SXK-THREE*SXK*SZK*SCT)
     >     +  SR6*(TWO*SZK-THREE*SZK*SZK*SCT))/SK2
      RETURN
      END
C*MODULE DFTXCC *DECK VS98XF
      SUBROUTINE VS98XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                  GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                  FTOTWT,XGRD,
     >                  VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA,
     >                  VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB,
     >                  TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IJZY)
C
C     REFERENCE
C     [A] T. V. VOORHIS & G. E. SCUSERIA, J. CHEM. PHYS. 109, 400 (1998)
C     [B] Y. ZHAO & D. G. TRUHLAR, J. CHEM. PHYS. 125, 194101 (2006)
C
      IMPLICIT NONE
C
C     CLOSED-SHELL OPEN-SHELL LOGICALS
C
      LOGICAL RHFGVB,UROHF
C
C     EXCHANGE ENERGY
C
      DOUBLE PRECISION XGRD
C
C     CHARGE DENSITY
C
      DOUBLE PRECISION ROA,ROB
C
C     CHARGE DENSITY GRADIENT AND COMPONENTS
C
      DOUBLE PRECISION GRDAA,GRDBB
      DOUBLE PRECISION GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB
C
C     KINETIC ENERGY DENSITY COMPONENTS
C
      DOUBLE PRECISION TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB
C
C     QUADRATURE WEIGHTS
C
      DOUBLE PRECISION FTOTWT
C
C     THE XC POTENTIAL
C
      DOUBLE PRECISION VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA
      DOUBLE PRECISION VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB
C
C     PARAMETER CONTROL
C
      INTEGER IJZY
C
C     WORKING VARIABLES
C
      DOUBLE PRECISION RRHO, RHO43, RHO13, RHOO, RHO53, RHO83
      DOUBLE PRECISION GAMMA
C
      DOUBLE PRECISION TAUN,TAUU
      DOUBLE PRECISION F43, F53, F83
      DOUBLE PRECISION GX, GG, X, Z, KX,XK,ZK
      DOUBLE PRECISION ONE, TWO
      DOUBLE PRECISION CF, AXLSDA, R1, R2, R3, R4, R5, R6
C
C     FUNCTIONAL DERIVATIVES
C
      DOUBLE PRECISION DXDR, DXDG, DZDR, DZDT, DGDX, DGDZ
C
      PARAMETER (CF  = 9.115599720D+00, AXLSDA = -0.9305257363491D+00)
      PARAMETER (GG  = 0.00186726D+00)
      PARAMETER (F43 = 4.0D+00/3.0D+00,F53=5.0D+00/3.0D+00)
      PARAMETER (F83 = 8.0D+00/3.0D+00)
      PARAMETER (ONE = 1.0D+00, TWO = 2.0D+00)
C
      IF (IJZY.EQ.1) THEN
C     PARAMETERS FOR VS98
        R1 = -9.800683D-01
        R2 = -3.556788D-03
        R3 =  6.250326D-03
        R4 = -2.354518D-05
        R5 = -1.282732D-04
        R6 =  3.574822D-04
      ELSEIF (IJZY.EQ.2) THEN
C     PARAMETERS FOR M06-L
        R1 =  6.012244D-01*AXLSDA
        R2 =  4.748822D-03*AXLSDA
        R3 = -8.635108D-03*AXLSDA
        R4 = -9.308062D-06*AXLSDA
        R5 =  4.482811D-05*AXLSDA
        R6 =  0.000000D+00
      ELSEIF (IJZY.EQ.3) THEN
C     PARAMETERS FOR M06-HF
        R1 = -1.179732D-01*AXLSDA
        R2 = -2.500000D-03*AXLSDA
        R3 = -1.180065D-02*AXLSDA
        R4 =  0.000000D+00
        R5 =  0.000000D+00
        R6 =  0.000000D+00
      ELSEIF (IJZY.EQ.4) THEN
C     PARAMETERS FOR M06
        R1 =  1.422057D-01*AXLSDA
        R2 =  7.370319D-04*AXLSDA
        R3 = -1.601373D-02*AXLSDA
        R4 =  0.000000D+00
        R5 =  0.000000D+00
        R6 =  0.000000D+00
      ENDIF
C
      IF (RHFGVB) THEN
        RHOO  = ROA
        RHO43 = RHOO**F43
        RRHO  = 1.0D+00/RHOO
        RHO13 = RHO43*RRHO
        RHO53 = RHOO**F53
        RHO83 = RHO53*RHOO
C
        TAUN = TAUXA + TAUYA + TAUZA
        TAUU = TAUN
C
        GAMMA = GRDAA
        X     = GAMMA/RHO83
        DXDR  =-F83*X*RRHO
        DXDG  = ONE/RHO83
        Z     = TAUU/RHO53 - CF
        DZDR  =-F53 * TAUU/RHO83
        DZDT  = ONE/RHO53
C
        KX = ONE + GG*X + GG*Z
        XK = X/KX
        ZK = Z/KX
C
        CALL GVT4(GX,DGDX,DGDZ,XK,ZK,KX,GG,GG,R1,R2,R3,R4,R5,R6)
C
        XGRD = XGRD + TWO*RHO43*GX*FTOTWT
C
C     FUNCTIONAL DERIVATIVES
C
        VXCA1 = VXCA1 + F43*RHO13*GX + RHO43*(DGDX*DXDR + DGDZ*DZDR)
C
        DUMAX = DUMAX + 2.0D+00*(RHO43*(DGDX*DXDG)*GRADXA)
        DUMAY = DUMAY + 2.0D+00*(RHO43*(DGDX*DXDG)*GRADYA)
        DUMAZ = DUMAZ + 2.0D+00*(RHO43*(DGDX*DXDG)*GRADZA)
C
        DMXCA = DMXCA + RHO43*(DGDZ*DZDT)
      END IF
C
      IF (UROHF) THEN
C
C     ALPHA
C
        RHOO  = ROA
        RHO43 = RHOO**F43
        RRHO  = 1.0D+00/RHOO
        RHO13 = RHO43*RRHO
        RHO53 = RHOO**F53
        RHO83 = RHO53*RHOO
C
        TAUN =(TAUXA + TAUYA + TAUZA)*(ONE/TWO)
        TAUU = TAUN*TWO
C
        GAMMA = GRDAA
C
        X    = GAMMA/RHO83
        DXDR = -F83*X*RRHO
        DXDG = ONE/RHO83
        Z    = TAUU/RHO53 - CF
        DZDR = -F53 * TAUU/RHO83
        DZDT = ONE/RHO53
C
        KX = ONE + GG*X + GG*Z
        XK = X/KX
        ZK = Z/KX
C
        CALL GVT4(GX,DGDX,DGDZ,XK,ZK,KX,GG,GG,R1,R2,R3,R4,R5,R6)
C
        XGRD = XGRD + RHO43*GX*FTOTWT
C
C     FUNCTIONAL DERIVATIVES
C
        VXCA1 = VXCA1 + F43*RHO13*GX +
     >                  RHO43*(DGDX*DXDR + DGDZ*DZDR)
C
        DUMAX = DUMAX + 2.0D+00*(RHO43*(DGDX*DXDG)*GRADXA)
        DUMAY = DUMAY + 2.0D+00*(RHO43*(DGDX*DXDG)*GRADYA)
        DUMAZ = DUMAZ + 2.0D+00*(RHO43*(DGDX*DXDG)*GRADZA)
C
        DMXCA = DMXCA + RHO43*(DGDZ*DZDT)
C
C     BETA
C
        RHOO  = ROB
        RHO43 = RHOO**F43
        RRHO  = 1.0D+00/RHOO
        RHO13 = RHO43*RRHO
        RHO53 = RHOO**F53
        RHO83 = RHO53*RHOO
C
        TAUN =(TAUXB + TAUYB + TAUZB)*(ONE/TWO)
        TAUU = TWO*TAUN
C
        GAMMA = GRDBB
C
        X    = GAMMA/RHO83
        DXDR = -F83*X*RRHO
        DXDG = ONE/RHO83
        Z    = TAUU/RHO53 - CF
        DZDR = -F53 * TAUU/RHO83
        DZDT = ONE/RHO53
C
        KX = ONE + GG*X + GG*Z
        XK = X/KX
        ZK = Z/KX
C
        CALL GVT4(GX,DGDX,DGDZ,XK,ZK,KX,GG,GG,R1,R2,R3,R4,R5,R6)
C
        XGRD = XGRD + RHO43*GX*FTOTWT
C
C     FUNCTIONAL DERIVATIVES
C
        VXCB1 = VXCB1 + F43*RHO13*GX +
     >                  RHO43*(DGDX*DXDR + DGDZ*DZDR)
C
        DUMBX = DUMBX + 2.0D+00*(RHO43*(DGDX*DXDG)*GRADXB)
        DUMBY = DUMBY + 2.0D+00*(RHO43*(DGDX*DXDG)*GRADYB)
        DUMBZ = DUMBZ + 2.0D+00*(RHO43*(DGDX*DXDG)*GRADZB)
C
        DMXCB = DMXCB + RHO43*(DGDZ*DZDT)
      ENDIF
      RETURN
      END
C*MODULE DFTXCC *DECK M06CF
      SUBROUTINE M06CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                 GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                 FTOTWT,ECF,
     >                 VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA,
     >                 VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB,
     >                 TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IJZY)
C
C     REFERENCES:
C     [A] ZHAO, Y. & TRUHLAR, D. G. J. CHEM. PHYS. 2006, 125, 194101;
C     [B] ZHAO, Y. & TRUHLAR, D. G. J. PHYS. CHEM. A (2006),110(49),
C         13126-13130.
C
      IMPLICIT NONE
C
C     CLOSED-SHELL OPEN-SHELL LOGICALS
C
      LOGICAL RHFGVB,UROHF
C
C     EXCHANGE ENERGY
C
      DOUBLE PRECISION ECF
C
C     CHARGE DENSITY
C
      DOUBLE PRECISION ROA,ROB
C
C     CHARGE DENSITY GRADIENT AND COMPONENTS
C
      DOUBLE PRECISION GRDAA,GRDBB
      DOUBLE PRECISION GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB
C
C     KINETIC ENERGY DENSITY COMPONENTS
C
      DOUBLE PRECISION TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB
C
C     QUADRATURE WEIGHTS
C
      DOUBLE PRECISION FTOTWT
C
C     THE XC POTENTIAL
C
      DOUBLE PRECISION VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA
      DOUBLE PRECISION VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB
C
C     PARAMETER CONTROL
C
      INTEGER IJZY
C
C     WORKING VARIABLES
C
      DOUBLE PRECISION F1, F2, F3, F4,COPP
      DATA F1/1.0D+00/,F2/2.0D+00/,F3/3.0D+00/,F4/4.0D+00/,
     >     COPP/0.0031D+00/
C
C     CALL TO THE M05CSS SUBROUTINE
C
      DOUBLE PRECISION PA,GAA,TA,FA,FPA,FGA,FTA,EUA,CHIA,EUPA,
     >                 CHIAP,CHIAG
      DOUBLE PRECISION PB,GBB,TB,FB,FPB,FGB,FTB,EUB,CHIB,EUPB,
     >                 CHIBP,CHIBG
C
      DOUBLE PRECISION  SOP, SOPP0, SOPP1,SOPP2, SOPP3, SOPP4
      DOUBLE PRECISION PI, PI34, F13,RS,RSP,ZETA,
     >                 DZDA,DZDB,POTLC,DLDS,DLDZ,D2LDSS,D2LDSZ,D2LDZZ
      DOUBLE PRECISION P, EUEG, U, W
      DOUBLE PRECISION DUDCHIA,DUDCHIB,DUDPA,DUDPB,DUDGA,DUDGB,
     >                 DWDU,DWDPA,DWDPB, DWDGA,DWDGB,EUEGPA,EUEGPB
C
      SOP=1.0D0
C
      IF (IJZY.EQ.1) THEN
C     PARAMETERS FOR M06-L CORRELATION
         SOPP0 =  6.042374D-01
         SOPP1 =  1.776783D+02
         SOPP2 = -2.513252D+02
         SOPP3 =  7.635173D+01
         SOPP4 = -1.255699D+01
      ELSEIF (IJZY.EQ.2) THEN
C     PARAMETERS FOR M06-HF CORRELATION
         SOPP0 =  1.674634D+00
         SOPP1 =  5.732017D+01
         SOPP2 =  5.955416D+01
         SOPP3 = -2.311007D+02
         SOPP4 =  1.255199D+02
      ELSEIF (IJZY.EQ.3) THEN
C     PARAMETERS FOR M06 CORRELATION
         SOPP0 =  3.741539D+00
         SOPP1 =  2.187098D+02
         SOPP2 = -4.531252D+02
         SOPP3 =  2.936479D+02
         SOPP4 = -6.287470D+01
      ELSEIF (IJZY.EQ.4) THEN
C     PARAMETERS FOR M06-2X CORRELATION
         SOPP0 =  8.833596D-01
         SOPP1 =  3.357972D+01
         SOPP2 = -7.043548D+01
         SOPP3 =  4.978271D+01
         SOPP4 = -1.852891D+01
      ENDIF
C
      CALL VS98CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >            GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >            FTOTWT,ECF,
     >            VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA,
     >            VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB,
     >            TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IJZY+1)
C
      PI   = F4*ATAN(F1)
      PI34 = F3 / (F4*PI)
      F13  = F1 / F3
C
      IF (RHFGVB) THEN
C
C     GET DENSITY, GRADIENT, & TAU FOR THE ALPHA SPIN FROM THE TOTAL
C
        PA  = ROA
        GAA = GRDAA
C
C     IN THE M05CSS SUBROUTINE, WE USE 2*TA AS THE TAU, SO WE DO NOT
C     DIVIDE THE TAU BY 2 HERE
C
        TA = TAUXA + TAUYA + TAUZA
C
        CALL M06CSS(PA,GAA,TA,FA,FPA,FGA,FTA,EUA,
     >              CHIA,EUPA,CHIAP,CHIAG,IJZY)
C
        PB    = PA
        GBB   = GAA
        TB    = TA
        FB    = FA
        FPB   = FPA
        FGB   = FGA
        FTB   = FTA
        EUB   = EUA
        CHIB  = CHIA
        EUPB  = EUPA
        CHIBP = CHIAP
        CHIBG = CHIAG
C
        ECF = ECF + 2.0D+00*FA*FTOTWT
C
        VXCA1 = VXCA1 + FPA
C
        DUMAX = DUMAX + 2.0D+00*(FGA*GRADXA)
        DUMAY = DUMAY + 2.0D+00*(FGA*GRADYA)
        DUMAZ = DUMAZ + 2.0D+00*(FGA*GRADZA)
C
        DMXCA = DMXCA + FTA
      END IF
      IF (UROHF) THEN
C
C       ALPHA
C
        PA  = ROA
        GAA = GRDAA
C
C  IN THE M05CSS SUBROUTINE, WE USE 2*TA AS THE TAU
C
        TA = TAUXA + TAUYA + TAUZA
C
        CALL M06CSS(PA,GAA,TA,FA,FPA,FGA,FTA,EUA,
     >              CHIA,EUPA,CHIAP,CHIAG,IJZY)
C
        ECF = ECF + FA*FTOTWT
C
        VXCA1 = VXCA1 + FPA
C
        DUMAX = DUMAX + 2.0D+00*(FGA*GRADXA)
        DUMAY = DUMAY + 2.0D+00*(FGA*GRADYA)
        DUMAZ = DUMAZ + 2.0D+00*(FGA*GRADZA)
C
        DMXCA = DMXCA + FTA
C
C       BETA
C
        PB  = ROB
        GBB = GRDBB
C
C     IN THE M05CSS SUBROUTINE, WE USE 2*TB AS THE TAU,
C
        TB = TAUXB + TAUYB + TAUZB
C
        CALL M06CSS(PB,GBB,TB,FB,FPB,FGB,FTB,EUB,
     >               CHIB,EUPB,CHIBP,CHIBG,IJZY)
C
        ECF = ECF + FB*FTOTWT
C
        VXCB1 = VXCB1 + FPB
C
        DUMBX = DUMBX + 2.0D+00*(FGB*GRADXB)
        DUMBY = DUMBY + 2.0D+00*(FGB*GRADYB)
        DUMBZ = DUMBZ + 2.0D+00*(FGB*GRADZB)
C
        DMXCB = DMXCB + FTB
      ENDIF
C
      P    = PA + PB
      RS   =(PI34/P) ** F13
      RSP  =-RS/(F3*P)
      ZETA =(PA-PB)/P
      DZDA =(F1-ZETA)/P
      DZDB =(-F1-ZETA)/P
C
      CALL LSDAC(RS,ZETA,POTLC,DLDS,DLDZ,D2LDSS,D2LDSZ,
     >           D2LDZZ)
C
      EUEG = P*POTLC - EUA - EUB
      U    = COPP*(CHIA+CHIB)/(F1 + COPP*(CHIA+CHIB))
      W    = SOPP0+U*(SOPP1+U*(SOPP2+U*(SOPP3+U*SOPP4)))
C
      ECF = ECF + SOP*EUEG*W*FTOTWT
C
      DUDCHIA = COPP/(F1 + COPP*(CHIA+CHIB))**2
      DUDCHIB = COPP/(F1 + COPP*(CHIA+CHIB))**2
C
      DUDPA = DUDCHIA*CHIAP
      DUDPB = DUDCHIB*CHIBP
      DUDGA = DUDCHIA*CHIAG
      DUDGB = DUDCHIB*CHIBG
      DWDU  = SOPP1+U*(F2*SOPP2+U*(F3*SOPP3+U*F4*SOPP4))
      DWDPA = DWDU*DUDPA
      DWDPB = DWDU*DUDPB
      DWDGA = DWDU*DUDGA
      DWDGB = DWDU*DUDGB
C
      EUEGPA = POTLC + P*DLDS*RSP + P*DLDZ*DZDA - EUPA
      EUEGPB = POTLC + P*DLDS*RSP + P*DLDZ*DZDB - EUPB
C
      IF (RHFGVB) THEN
        VXCA1 = VXCA1 + SOP*(EUEGPA*W + EUEG*DWDPA)
C
        DUMAX = DUMAX + 2.0D+00*SOP*EUEG*DWDGA*GRADXA
        DUMAY = DUMAY + 2.0D+00*SOP*EUEG*DWDGA*GRADYA
        DUMAZ = DUMAZ + 2.0D+00*SOP*EUEG*DWDGA*GRADZA
      END IF
C
      IF (UROHF) THEN
        VXCA1 = VXCA1 + SOP*(EUEGPA*W + EUEG*DWDPA)
C
        VXCB1 = VXCB1 + SOP*(EUEGPB*W + EUEG*DWDPB)
C
        DUMAX = DUMAX + 2.0D+00*SOP*EUEG*DWDGA*GRADXA
        DUMAY = DUMAY + 2.0D+00*SOP*EUEG*DWDGA*GRADYA
        DUMAZ = DUMAZ + 2.0D+00*SOP*EUEG*DWDGA*GRADZA
C
        DUMBX = DUMBX + 2.0D+00*SOP*EUEG*DWDGB*GRADXB
        DUMBY = DUMBY + 2.0D+00*SOP*EUEG*DWDGB*GRADYB
        DUMBZ = DUMBZ + 2.0D+00*SOP*EUEG*DWDGB*GRADZB
      ENDIF
      RETURN
      END
C
C*MODULE DFTXCC *DECK M06CSS
      SUBROUTINE M06CSS(PX,GX,TX,F,FP,FG,FT,EUEG,CHI,EUEGP,
     >                  CHIP,CHIG,IJZY)
      IMPLICIT NONE
C
C     COMPUTE THE SAME-SPIN PART OF THE M05 CORRELATION FUNCTIONAL
C     FOR ONE GRID POINT AND ONE SPIN-CASE.
C
      INTEGER IJZY
      DOUBLE PRECISION PX, GX, TX, F, FP, FG, FT
      DOUBLE PRECISION EUEG, CHI, EUEGP, CHIP, CHIG
      DOUBLE PRECISION PT25, F1, F2, F3, F4, F8
      DOUBLE PRECISION SS, SSS0,SSS1, SSS2, SSS3, SSS4, CSS
      DOUBLE PRECISION PI, PI34, F13, F83
      DOUBLE PRECISION RS, D, FSCC, RSP, DFSCCP, DFSCCG
      DOUBLE PRECISION E, W, U, DFSCCT, DUDCHI, DWDU, DWDP, DWDG
      DOUBLE PRECISION D2LDSS,D2LDSZ,D2LDZZ,POTLC,DLDS,DLDZ
C
      DATA PT25/0.25D+00/, F1/1.0D+00/, F2/2.0D+00/, F3/3.0D+00/,
     >     CSS/0.06D+00/ , F4/4.0D+00/, F8/8.0D+00/
C
      SS=1.0D+00
C
      IF (IJZY.EQ.1) THEN
C     PARAMETERS FOR M06-L CORRELATION
         SSS0 =  5.349466D-01
         SSS1 =  5.396620D-01
         SSS2 = -3.161217D+01
         SSS3 =  5.149592D+01
         SSS4 = -2.919613D+01
      ELSEIF (IJZY.EQ.2) THEN
C     PARAMETERS FOR M06-HF CORRELATION
         SSS0 =  1.023254D-01
         SSS1 = -2.453783D+00
         SSS2 =  2.913180D+01
         SSS3 = -3.494358D+01
         SSS4 =  2.315955D+01
      ELSEIF (IJZY.EQ.3) THEN
C     PARAMETERS FOR M06 CORRELATION
         SSS0 =  5.094055D-01
         SSS1 = -1.491085D+00
         SSS2 =  1.723922D+01
         SSS3 = -3.859018D+01
         SSS4 =  2.845044D+01
      ELSEIF (IJZY.EQ.4) THEN
C     PARAMETERS FOR M06-2X CORRELATION
         SSS0 =  3.097855D-01
         SSS1 = -5.528642D+00
         SSS2 =  1.347420D+01
         SSS3 = -3.213623D+01
         SSS4 =  2.846742D+01
      ENDIF
C
      PI    = F4*ATAN(F1)
      PI34  = F3 / (F4*PI)
      F13   = F1 / F3
      F83   = F8 / F3
      RS    =(PI34/PX) ** F13
C
      CALL LSDAC(RS,F1,POTLC,DLDS,DLDZ,D2LDSS,D2LDSZ,D2LDZZ)
C
      EUEG = PX*POTLC
      D    = TX - PT25*GX/PX
C
      CHI  = GX/PX**F83
      U    = CSS*CHI/(F1 + CSS*CHI)
      W    = SSS0+U*(SSS1+U*(SSS2+U*(SSS3+U*SSS4)))
      FSCC = D/TX
      E    = FSCC*W*EUEG
      F    = E*SS
      RSP  =-RS/(F3*PX)
      CHIG = F1/PX**F83
      CHIP =-F83*CHI/PX
C
      DFSCCP = PT25*GX/(TX*PX**2)
      DFSCCG =-PT25/(TX*PX)
      DFSCCT = PT25*GX/(PX*TX**2)
      DUDCHI = CSS/((F1+CSS*CHI)**2)
      DWDU   = SSS1+U*(F2*SSS2+U*(F3*SSS3+U*F4*SSS4))
      DWDP   = DWDU*DUDCHI*CHIP
      DWDG   = DWDU*DUDCHI*CHIG
C
      EUEGP = POTLC + PX*DLDS*RSP
C
      FP = SS*(DFSCCP*W*EUEG
     >   + FSCC*DWDP*EUEG
     >   + FSCC*W*EUEGP)
      FG = SS*(DFSCCG*W*EUEG
     >   + FSCC*DWDG*EUEG)
      FT = SS*(DFSCCT*W*EUEG)
      RETURN
      END
C*MODULE DFTXCC *DECK M06XF
      SUBROUTINE M06XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                 GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                 FTOTWT,XGRD,
     >                 VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA,
     >                 VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB,
     >                 TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IJZY)
C
C     REFERENCES:
C     [A] ZHAO, Y. & TRUHLAR, D. G. J. CHEM. PHYS. 2006, 125, 194101;
C     [B] ZHAO, Y. & TRUHLAR, D. G. J. PHYS. CHEM. A (2006),110(49),
C         13126-13130.
C
      IMPLICIT NONE
C
C     CLOSED-SHELL OPEN-SHELL LOGICALS
C
      LOGICAL RHFGVB,UROHF
C
C     EXCHANGE ENERGY
C
      DOUBLE PRECISION XGRD
C
C     CHARGE DENSITY
C
      DOUBLE PRECISION ROA,ROB
C
C     CHARGE DENSITY GRADIENT AND COMPONENTS
C
      DOUBLE PRECISION GRDAA,GRDBB
      DOUBLE PRECISION GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB
C
C     KINETIC ENERGY DENSITY COMPONENTS
C
      DOUBLE PRECISION TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB
C
C     QUADRATURE WEIGHTS
C
      DOUBLE PRECISION FTOTWT
C
C     THE XC POTENTIAL
C
      DOUBLE PRECISION VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA
      DOUBLE PRECISION VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB
C
C     PARAMETER CONTROL
C
      INTEGER IJZY
C
C     WORKING VARIABLES
C
      DOUBLE PRECISION AT1, AT2, AT3, AT4, AT5, AT6, AT7, AT8, AT9
      DOUBLE PRECISION AT, AT10, AT11, AT0, C1, C2, FL, FNL
      DOUBLE PRECISION RRHO, RHO43, RHO13, RHOO, RHO53
      DOUBLE PRECISION GAMMA2, GAMMA
      DOUBLE PRECISION TAUUEG, TSIG, WSIG, W1, W2, W3, W4, W5, W6
      DOUBLE PRECISION W7, W8, W9, W10, W11, FSIG, PI
C
      DOUBLE PRECISION TAUN,TAUU
      DOUBLE PRECISION F23, F53, F13
      DOUBLE PRECISION F4O3
      DOUBLE PRECISION ONE, TWO, THREE, FOUR, FIVE, SIX, SEVEN, EIGHT
      DOUBLE PRECISION NINE, F10, F11
      DOUBLE PRECISION CS, AX, P32, S, S2, EN, ED, E, DE, DEN, DED
C
C     FUNCTIONAL DERIVATIVES
C
      DOUBLE PRECISION DFDW, DWDT, DTDR, DTDTAU, DGGADR, DFDR
      DOUBLE PRECISION DFDTAU, DGGADG
C
      PARAMETER(PI=3.1415926535897932384626433832795D+00)
C
      PARAMETER(F13=1.0D+00/3.0D+00,F4O3=4.0D+00/3.0D+00)
      PARAMETER(F23=2.0D+00/3.0D+00, F53=5.0D+00/3.0D+00)
C
      PARAMETER( ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00,  FOUR=4.0D+00,
     >          FIVE=5.0D+00, SIX=6.0D+00, SEVEN=7.0D+00, EIGHT=8.0D+00,
     >          NINE=9.0D+00, F10=10.D+00,   F11=11.D+00)
C
      IF (IJZY.EQ.1) THEN
C     PARAMETERS FOR M06-L
        AT0  =  3.987756D-01
        AT1  =  2.548219D-01
        AT2  =  3.923994D-01
        AT3  = -2.103655D+00
        AT4  = -6.302147D+00
        AT5  =  1.097615D+01
        AT6  =  3.097273D+01
        AT7  = -2.318489D+01
        AT8  = -5.673480D+01
        AT9  =  2.160364D+01
        AT10 =  3.421814D+01
        AT11 = -9.049762D+00
       ELSEIF (IJZY.EQ.2) THEN
C     PARAMETERS FOR M06-HF
        AT0  =  1.179732D-01
        AT1  = -1.066708D+00
        AT2  = -1.462405D-01
        AT3  =  7.481848D+00
        AT4  =  3.776679D+00
        AT5  = -4.436118D+01
        AT6  = -1.830962D+01
        AT7  =  1.003903D+02
        AT8  =  3.864360D+01
        AT9  = -9.806018D+01
        AT10 = -2.557716D+01
        AT11 =  3.590404D+01
       ELSEIF (IJZY.EQ.3) THEN
C     PARAMETERS FOR M06
        AT0  =  5.877943D-01
        AT1  = -1.371776D-01
        AT2  =  2.682367D-01
        AT3  = -2.515898D+00
        AT4  = -2.978892D+00
        AT5  =  8.710679D+00
        AT6  =  1.688195D+01
        AT7  = -4.489724D+00
        AT8  = -3.299983D+01
        AT9  = -1.449050D+01
        AT10 =  2.043747D+01
        AT11 =  1.256504D+01
       ELSEIF (IJZY.EQ.4) THEN
C     PARAMETERS FOR M06-2X
        AT0  =  4.600000D-01
        AT1  = -2.206052D-01
        AT2  = -9.431788D-02
        AT3  =  2.164494D+00
        AT4  = -2.556466D+00
        AT5  = -1.422133D+01
        AT6  =  1.555044D+01
        AT7  =  3.598078D+01
        AT8  = -2.722754D+01
        AT9  = -3.924093D+01
        AT10 =  1.522808D+01
        AT11 =  1.522227D+01
      ENDIF
C
      AT = 1.0D+00
      IF (IJZY.LT.4) CALL VS98XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                    GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                    FTOTWT,XGRD,
     >                    VXCA1,DUMAX,DUMAY,DUMAZ,DMXCA,
     >                    VXCB1,DUMBX,DUMBY,DUMBZ,DMXCB,
     >                    TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IJZY+1)
C
      C1 = 0.2195149727645171D+00
      C2 = C1/0.804D+00
C
C     SCALE FACTORS FOR LOCAL AND NON-LOCAL CONTRIBUTIONS.
C
      FL  = 1.0D+00
      FNL = 1.0D+00
      CS  = 0.5D+00/(3.0D+00*PI*PI)**F13
      P32 =(3.D+00*PI**2)**F23
      AX  =(-0.75D+00)*(3.0D+00/PI)**F13
C
      IF (RHFGVB)THEN
        RHOO  = TWO*ROA
        RHO43 = RHOO**F4O3
        RRHO  = 1D+00/RHOO
        RHO13 = RHO43*RRHO
        RHO53 = RHOO**F53
C
        TAUN   = TAUXA + TAUYA + TAUZA
        TAUU   = TAUN
        TAUUEG = 0.3D+00*P32*RHO53
        TSIG   = TAUUEG/TAUU
        WSIG   =(TSIG-ONE)/(TSIG+ONE)
        W1     = WSIG
        W2     = WSIG*W1
        W3     = WSIG*W2
        W4     = WSIG*W3
        W5     = WSIG*W4
        W6     = WSIG*W5
        W7     = WSIG*W6
        W8     = WSIG*W7
        W9     = WSIG*W8
        W10    = WSIG*W9
        W11    = WSIG*W10
        FSIG   = AT*(AT0 + AT1*W1 + AT2*W2 + AT3*W3
     >                   + AT4*W4 + AT5*W5 + AT6*W6   + AT7*W7
     >                   + AT8*W8 + AT9*W9 + AT10*W10 + AT11*W11)
C
        GAMMA2 = 4.0D+00*GRDAA
        GAMMA  = SQRT(GAMMA2)
        S      = CS*GAMMA/RHO43
        S2     = S*S
        EN     = C1*S2
        ED     = ONE + C2*S2
        E      = EN/ED
C
        XGRD = XGRD + RHO43*AX*(FL+FNL*E)*FSIG*FTOTWT
C
C     FUNCTIONAL DERIVATIVES
C
        DEN    = TWO*C1*S
        DED    = TWO*C2*S
        DE     = (DEN*ED-EN*DED)/(ED*ED)
        DFDW   = AT*(AT1 + TWO*AT2*W1 + THREE*AT3*W2
     >         + FOUR*AT4*W3 + FIVE*AT5*W4
     >         + SIX*AT6*W5 + SEVEN*AT7*W6
     >         + EIGHT*AT8*W7 + NINE*AT9*W8
     >         + F10*AT10*W9 + F11*AT11*W10)
        DWDT   = TWO/((ONE + TSIG)**2)
        DTDR   = (0.5D+00*P32*RHO13*RHO13)/TAUU
        DTDTAU =-TAUUEG/TAUU**2
        DGGADR = F4O3*RHO13*AX*(FL+FNL*(E-S*DE))
        DFDR   = DFDW*DWDT*DTDR
        DFDTAU = DFDW*DWDT*DTDTAU
        DGGADG =(FNL*DE*S/(TWO*GAMMA2))
C
        VXCA1 = VXCA1 + DGGADR*FSIG+(FL+FNL*E)*AX*RHO43*DFDR
C
        DUMAX = DUMAX + 2.0D+00*TWO*DGGADG*AX*RHO43*FSIG*GRADXA
        DUMAY = DUMAY + 2.0D+00*TWO*DGGADG*AX*RHO43*FSIG*GRADYA
        DUMAZ = DUMAZ + 2.0D+00*TWO*DGGADG*AX*RHO43*FSIG*GRADZA
C
        DMXCA = DMXCA + 0.5D+00*RHO43*AX*(FL+FNL*E)*DFDTAU
      END IF
C
      IF (UROHF) THEN
C
C  USE SPIN DENSITY FUNCTIONAL THEORY IE N-->2N
C                               EX=(1/2)EX[2*ALPHA] + (1/2)EX[2*BETA]
C
C     ALPHA
C
        RHOO  = TWO*ROA
        RHO43 = RHOO**F4O3
        RRHO  = 1.0D+00/RHOO
        RHO13 = RHO43*RRHO
        RHO53 = RHOO**F53
C
        TAUN   = (TAUXA + TAUYA + TAUZA)*(ONE/TWO)
        TAUU   = TWO*TAUN
        TAUUEG = 0.3D+00*P32*RHO53
        TSIG   = TAUUEG/TAUU
        WSIG   =(TSIG-ONE)/(TSIG+ONE)
        W1     = WSIG
        W2     = WSIG*W1
        W3     = WSIG*W2
        W4     = WSIG*W3
        W5     = WSIG*W4
        W6     = WSIG*W5
        W7     = WSIG*W6
        W8     = WSIG*W7
        W9     = WSIG*W8
        W10    = WSIG*W9
        W11    = WSIG*W10
C
        FSIG = AT*(AT0    + AT1*W1 + AT2*W2   + AT3*W3
     >       +     AT4*W4 + AT5*W5 + AT6*W6   + AT7*W7
     >       +     AT8*W8 + AT9*W9 + AT10*W10 + AT11*W11)
C
        GAMMA2 = GRDAA
        GAMMA2 = FOUR*GAMMA2
        GAMMA  = SQRT(GAMMA2)
C
        S      = CS*GAMMA/RHO43
        S2     = S*S
        EN     = C1*S2
        ED     = ONE + C2*S2
        E      = EN/ED
C
        XGRD = XGRD + RHO43*AX*(FL+FNL*E)*FSIG*0.5D+00*FTOTWT
C
C     FUNCTIONAL DERIVATIVES
C
        DEN   = TWO*C1*S
        DED   = TWO*C2*S
        DE    =(DEN*ED-EN*DED)/(ED*ED)
C
        DFDW = AT*(AT1 + TWO*AT2*W1 + THREE*AT3*W2
     >       +     FOUR*AT4*W3      + FIVE*AT5*W4
     >       +     SIX*AT6*W5       + SEVEN*AT7*W6
     >       +     EIGHT*AT8*W7     + NINE*AT9*W8
     >       +     F10*AT10*W9      + F11*AT11*W10)
C
        DWDT   = TWO/((ONE + TSIG)**2)
        DTDR   = TWO*(0.5D+00*P32*RHO13*RHO13)/TAUU
        DTDTAU =-TWO*TAUUEG/TAUU**2
        DGGADR = TWO*F4O3*RHO13*AX*(FL+FNL*(E-S*DE))
        DFDR   = DFDW*DWDT*DTDR
        DFDTAU = DFDW*DWDT*DTDTAU
        DGGADG = FOUR*(FNL*DE*S/(TWO*GAMMA2))
C
        VXCA1 = VXCA1 + (DGGADR*FSIG+(FL+FNL*E)*AX*RHO43*DFDR)*0.5D+00
C
        DUMAX = DUMAX + 2.0D+00*DGGADG*AX*RHO43*FSIG*0.5D+00*GRADXA
        DUMAY = DUMAY + 2.0D+00*DGGADG*AX*RHO43*FSIG*0.5D+00*GRADYA
        DUMAZ = DUMAZ + 2.0D+00*DGGADG*AX*RHO43*FSIG*0.5D+00*GRADZA
C
        DMXCA = DMXCA + 0.5D+00*0.5D+00*RHO43*AX*(FL+FNL*E)*DFDTAU
C
C     BETA
C
        RHOO  = TWO*ROB
        RHO43 = RHOO**F4O3
        RRHO  = 1.0D+00/RHOO
        RHO13 = RHO43*RRHO
        RHO53 = RHOO**F53
C
        TAUN   = (TAUXB + TAUYB + TAUZB)*(ONE/TWO)
        TAUU   = TWO*TAUN
        TAUUEG = 0.3D+00*P32*RHO53
        TSIG   = TAUUEG/TAUU
        WSIG   =(TSIG-ONE)/(TSIG+ONE)
        W1     = WSIG
        W2     = WSIG*W1
        W3     = WSIG*W2
        W4     = WSIG*W3
        W5     = WSIG*W4
        W6     = WSIG*W5
        W7     = WSIG*W6
        W8     = WSIG*W7
        W9     = WSIG*W8
        W10    = WSIG*W9
        W11    = WSIG*W10
C
        FSIG =AT*(AT0    + AT1*W1 + AT2*W2   + AT3*W3
     >       +    AT4*W4 + AT5*W5 + AT6*W6   + AT7*W7
     >       +    AT8*W8 + AT9*W9 + AT10*W10 + AT11*W11)
C
        GAMMA2 = GRDBB
        GAMMA2 = FOUR*GAMMA2
        GAMMA  = SQRT(GAMMA2)
        S      = CS*GAMMA/RHO43
        S2     = S*S
        EN     = C1*S2
        ED     = ONE + C2*S2
        E      = EN/ED
C
        XGRD = XGRD + RHO43*AX*(FL+FNL*E)*FSIG*0.5D+00*FTOTWT
C
C     FUNCTIONAL DERIVATIVES
C
        DEN   = TWO*C1*S
        DED   = TWO*C2*S
        DE    =(DEN*ED-EN*DED)/(ED*ED)
C
        DFDW = AT*(AT1 + TWO*AT2*W1 + THREE*AT3*W2
     >       +     FOUR*AT4*W3      + FIVE*AT5*W4
     >       +     SIX*AT6*W5       + SEVEN*AT7*W6
     >       +     EIGHT*AT8*W7     + NINE*AT9*W8
     >       +     F10*AT10*W9      + F11*AT11*W10)
C
        DWDT   = TWO/((ONE + TSIG)**2)
        DTDR   = TWO*(0.5D+00*P32*RHO13*RHO13)/TAUU
        DTDTAU =-TWO*TAUUEG/TAUU**2
        DGGADR = TWO*F4O3*RHO13*AX*(FL+FNL*(E-S*DE))
        DFDR   = DFDW*DWDT*DTDR
        DFDTAU = DFDW*DWDT*DTDTAU
        DGGADG = FOUR*(FNL*DE*S/(TWO*GAMMA2))
C
        VXCB1 = VXCB1 + (DGGADR*FSIG + (FL+FNL*E)*AX*RHO43*DFDR)*0.5D+00
C
        DUMBX = DUMBX + 2.0D+00* DGGADG*AX*RHO43*FSIG*0.5D+00*GRADXB
        DUMBY = DUMBY + 2.0D+00* DGGADG*AX*RHO43*FSIG*0.5D+00*GRADYB
        DUMBZ = DUMBZ + 2.0D+00* DGGADG*AX*RHO43*FSIG*0.5D+00*GRADZB
C
        DMXCB = DMXCB + 0.5D+00*0.5D+00*RHO43*AX*(FL+FNL*E)*DFDTAU
      ENDIF
      RETURN
      END

