C 11 APR 08 - HPTI - INPGDFT: PRINT CHANGES FOR NEW ANGULAR GRIDS
C  4 MAR 08 - SOK - META-GGA FUNCTIONALS (NAME CHANGE FROM DFTEXC.SRC)
C 20 AUG 07 - MC  - FIX A BUG IN LC CASE
C 24 MAR 07 - SS  - PW91,PBE,PBE0,X3LYP,B3LYP1,VWN1,OPTX FUNCTIONALS
C 24 MAR 07 - RP  - ADD P86 AND PZ81 CORRELATION FUNCTIONALS
C 22 DEC 06 - ST,NK,MC - ADD LC EXCHANGE FUNCTIONALS
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 19 SEP 05 - MWS - ADD TRUE NUCLEAR CHARGE ARRAY TO INFOA COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 16 JAN 04 - MWS - SYNCH COMMON BLOCK DFGRID
C 24 JAN 02 - MWS - FIX CHAR VARIABLE USAGE, CHANGE CORR.ONLY KEYWORDS
C                   TO VWN/LYP/OP, FIX LYP FUNCTIONAL FOR NB=0 CASE
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 13 JUN 01 - TT,SY,MK,DGF IMPLEMENT GRID-BASED DFT
C
C***********************************************************************
C     HOW TO ADD A NEW FUNCTIONAL INTO THE DFT GRID MODULE: 03/18/2007
C     BY SAROM SOK (AKA MR DFT), DEPT. OF CHEM., IOWA STATE UNIVERSITY.
C
C     STEP 1) MODIFY THE INPGDFT ROUTINE
C
C     YOU WILL NEED TO BUMP UP THE PARAMTER -NEX- OR -NCOR-
C     DEPENDING ON IF YOU ARE ADDING AN EXCHANGE OR CORRELATION
C     FUNCTIONAL.  THE ADDITION OF HYBRID FUNCTIONALS WILL NOT NEED TO
C     MODIFY THIS PARAMETER.
C
C     NOTICE THAT THE FOLLOWING ARRAYS ARE DEPENDING ON THE ABOVE
C     PARAMETER CHANGES:
C
C           -EX-, -COR- -EXNM-, -CORNM-, -NUMEX-, -NUMCOR-, -DOP-
C
C     CREATE A SPOT FOR YOUR FUNCTIONAL IN THE -EX- OR -COR- NAME ARRAY.
C     NOTE THAT IF YOU ADDING IN A HYBRID FUNCTIONAL AND ALL THE
C     INGREDIENTS OF THE HYBRID CURRENTLY EXISTS THEN YOU DO NOT NEED TO
C     MODIFY THE -EX- OR -COR- NAME ARRAY.
C     MODIFYING THE -EX- OR -COR- NAME ARRAYS ARE FOR NEW FUNCTIONALS
C
C     MODIFYING THE -EX- OR -COR- NAME ARRAY ALLOWS FOR THE MIXING AND
C     MATCHING OF EXCHANGE AND CORRELATION FUNCTIONALS IN DFTTYP.
C     FOR INSTANCE: DFTTYP = [-EX- NAME][-COR- NAME] = OLYP
C     WILL MATCH OPTX EXCHANGE WITH LYP CORRELATION.
C
C     ANOTHER EXAMPLE IS: DFTTYP = PW91PW91
C     WILL MATCH PW91 EXCHANGE WITH PW91 CORRELATION.  BUT THIS IS
C     ALREAD DEFINED AS DFTTYP = PW91
C
C     ONE MORE EXAMPLE: DFTTYP = BPW91
C     WHILL MATCH B88 EXCHANGE WITH PW91 CORRELATION
C
C     IF YOUR ADDING IN A NEW SET OF FUNCTIONALS (EXCHANGE+CORRELATION)
C     PLEASE MAKE SURE THAT THE INDEX OF THE EXCHANGE FUNCTIONAL IN
C     THE -EX- NAME ARRAY IS ALSO THE SAME AS THE INDEX OF THE
C     CORRELATION FUNCTIONAL IN -COR- NAME ARRAY.
C
C     EG: PW91 - EX[3] = PW91 AND COR[3] = PW91
C
C     THE -NEXFG- AND -NCORFG- CONTROL VALUES SHOULD ALSO CORRESPOND TO
C     THE INDEX OF THE FUNCTIONAL IN THE NAME ARRAYS.
C
C     EG: PW91 EXCHANGE: NEXFG = 3, EX[3] = PW91
C     EG: PW91 CORRELATION: NCORFG = 3, COR[3] = PW91
C
C     HYBRID EXCHANGE FUCTIONALS AND HYBRID CORRELATION FUNCTIONALS ARE
C     TO USE NEGATIVE CONTROL NUMBERS FOR -NEXFG- AND -NCORFG-.
C
C     ADD A NEW DATA ENTRY FOR YOUR FUNCTIONAL
C     THIS DATA ENTRY CORRESPONDS TO WHAT THE USER NEEDS TO PROVIDE
C     IN DFTTYP.
C
C     IF YOU ARE ADDING AN EXCHANGE FUNCTIONAL THEN YOU WILL ALSO NEED
C     TO MAKE AN ENTRY FOR YOUR FUNCTIONAL IN DATA DQOP. THE VALUE YOU
C     WILL PROVIDE IS ZERO.
C
C     NEW HYBRID FUNCTIONALS WILL NEED A COMMON BLOCK TO STORE THE
C     FACTORS FOR EACH OF ITS FUNCTIONALS. YOU ARE PERMITTED TO CREATE A
C     NEW COMMON BLOCK FOR YOUR HYBRID.
C
C     EG: COMMON /CX3LYP/ CPW91
C     EG: COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
C
C     ADD AN ENTRY FOR YOUR FUNCTIONALS IN THE SERIES OF IF STATEMENTS.
C
C     EXAMPLE:
C      ELSEIF(DFTYPE.EQ.OPTX)  THEN
C         NEXFG     = 7
C         EXCHR     =  EXNM(7)      <-- NAME ENTRY OF EXCHANGE NAME
C         CORCHR    = 'NONE'        <-- MANUAL ENTRY OF CORRELATION NAME
C         CSLT      = 1.05151D+00   <-- FACTORS FOR EACH FUNCTIONAL
C         COPTX     = 1.43169D+00
C
C     THE AMOUNT OF HATREE-FOCK EXCHANGE IS DEFINED IN -DFTTYP(3)-
C
C     NOTE THAT THE B88 EXCHANGE FUNCTIONAL ALREADY HAS SLATER-DIRAC
C     EXCHANGE (ALSO REFERRED TO HAS LDA EXCHANGE).
C
C     SO BE CAREFUL WHEN DEFINING -CSLT- AND -CB88- AT THE SAME TIME!
C
C     STEP 2) MODIFY THE CALCEXC ROUTINE
C
C     ADD AN ELSEIF ENTRY FOR YOUR FUNCTIONAL.
C
C     STEP 3) ADD YOUR FUNCTIONAL ROUTINE
C
C     A WORD OF CAUTION IF YOU DECIDE TO TAKE A STEAL A PARAMETER LIST
C     FROM AN EXISTING SUBROUTINE.
C
C     IF YOU COMPARE THE PARAMETER NAMES THE OF B88 EXCHANGE IN THE
C     ROUTINE DEFINITION TO THE PARAMETER NAMES IN THE ROUTINE CALL TO
C     THE B88 EXCHANGE IN CALCEXC THEN YOU WILL NOTICED THAT THEY ARE
C     DIFFERENT.
C
C***********************************************************************
C     SPECIAL COMMENTS
C
C     IN ADDING THE OPTX EXCHANGE FUNCTIONAL INTO GAMESS IT WAS
C     DISCOVERED THAT ADDITIONAL PREFACTOR (AKA COEEFICIENTS)
C     ASSIGNMENTS HAD TO BE DECLARED IN THE ROUTINE NAMEXC.
C
C     THE -CSLT- PREFACTOR IS USED BY OTHER FUNCTIONALS AND HAS AN
C     INITIAL VALUE OF ONE IN THE ROUTINE INPGDFT.  IN THE OPTX
C     EXCHANGE FUNCTIONAL THIS VALUE IS 1.05151D+00.  NOW IN ORDER TO
C     ALLOW THE OPTX FUNCTIONAL TO BE MIXED WITH THE AVAILABLE
C     CORRELATION FUNCTIONALS IN GAMESS THE -CSLT- PREFACTOR HAD TO BE
C     ASSIGNED IN THE NAMEXC ROUTINE.  IT IS IN NAMEXC THAT THE NEXFG-
C     AND -NCORFG- VALUES ARE DETERMINED FOR AD-HOC EXCHANGE/CORRELATION
C     FUNCTIONAL COMBINATIONS.
C
C     02/21/2008
C     META-GGA CAPABILITIES HAVE BEEN IMPLEMENTED INTO THE DFT GRID
C     MODULE.  THE FOLLOWING FUNCTIONALS ARE META-GGAS (TAU-DEPENDENT)
C       - M05, M05-2X, VS98, M06, M06-L, M06-2X, M06-HF -
C
C     BOTH ENERGIES AND GRADIENTS WITH META-GGAS ARE AVAILABLE.
C     HESSIANS ARE NOT YET AVAILABLE.
C***********************************************************************
C
C*MODULE DFTXCA  *DECK B88XF
      SUBROUTINE B88XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                   GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                   FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >                   XALF,B88X,DUMA,DUMXA,DUMYA,DUMZA,
     >                   DUMB,DUMXB,DUMYB,DUMZB)
C***********************************************************************
C     BECKE 88 EXCHANGE FUNCTIONAL
C
C     A. D. BECKE
C     PHYS. REV. A 38,3098 (1988)
C
C     -----
C     INPUT
C     -----
C
C     ROA,B   .... SPINDENSITIES
C     GRAD**  .... DENSITY GRADIENTS
C     TOTWT   .... VOLUME OF THIS GRID
C
C     ------
C     OUTPUT
C     ------
C
C     XALF    ....  EXCHANGE ENERGY OF LDA PART
C     B88X    ....  EXCHANGE ENERGY OF GGA PART
C     DUM***  ....  FIRST DERIVATIVES OF THE EXCHANGE FUNCTINALS
C     SIGKA,B ....  EXC=-1/2*INT(SIGK*RHO^(4/3))
C     DKDXA,B ....  DSIGK/DX
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00)
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      LOGICAL LCFLAG,LRINT
C***********************************************************************
       X1 = (3.0D+00/PI)**(ONE/3.0D+00)
       BECKEB=0.0042D+00
       FACE  =-1.5D+00*X1*(ONE/TWO )**(TWO/THREE)
       FACP  = 4.0D+00/3.0D+00
C***********************************************************************
C      CLOSED SHELL CALCULATION
C***********************************************************************
       IF (RHFGVB) THEN
         RHO=ROA
         GRD=GRDAA
         RHO13=RHO**(ONE/3.0D+00)
         RHOE =RHO13**4
         X    =SQRT(GRD)/RHOE
         Q    =SQRT(ONE+X*X)
         S    =LOG(X+Q)
         T    =ONE+6.0D+00*BECKEB*X*S
         G    =-BECKEB*X*X/T
         H    =(-TWO*BECKEB*X+6.0D+00*(BECKEB*X)**2*(-S+X/Q))/T**2
         R    =FACP*RHO13
         E    =FACE
         XCL  =E * RHOE
         XCNL =G * RHOE
         XALFA =CB88*TWO*XCL *FTOTWT
         B88XA =CB88*TWO*XCNL*FTOTWT
         U     =CB88*R*(E+G-X*H)
         PDUMA =U
         W     =CB88*0.5D+00*H/SQRT(GRD)
         PDUMXA=TWO*W*GRADXA
         PDUMYA=TWO*W*GRADYA
         PDUMZA=TWO*W*GRADZA
C     ----- FOR OP CORRELATION -----
         SIGKA= -(FACE+G)*TWO
         DKDXA= -H*TWO
C     ------------------------------
         IF (LCFLAG) THEN
            CALL LRCFUN(UROHF,ROA,ROB,GRDAA,GRDBB,
     >           GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >           FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >           XALFA,XALFB,B88XA,B88XB,PDUMA,PDUMXA,PDUMYA,
     >           PDUMZA,PDUMB,PDUMXB,PDUMYB,PDUMZB,FACE)
         END IF
         XALF   = XALF + XALFA
         B88X   = B88X + B88XA
         DUMA   = DUMA + PDUMA
         DUMXA  = DUMXA+ PDUMXA
         DUMYA  = DUMYA+ PDUMYA
         DUMZA  = DUMZA+ PDUMZA
C***********************************************************************
C      OPEN SHELL CALCULATION
C***********************************************************************
       ELSE IF (UROHF) THEN
C     ----- ALPHA CONTRIBTUTION ----
         RHO13A=ROA**(ONE/3.0D+00)
         RHOEA =RHO13A**4
         XA    =SQRT(GRDAA)/RHOEA
         QA    =SQRT(ONE+XA**2)
         SA    =LOG(XA+QA)
         TA    =ONE+6.0D+00*BECKEB*XA*SA
         GA    =-BECKEB*XA*XA/TA
         HA    =(-TWO*BECKEB*XA
     >            +6.0D+00*(BECKEB*XA)**2*(-SA+XA/QA))/TA**2
         RA    =FACP*RHO13A
         EA    =FACE
         UA     =CB88*RA*(EA+GA-XA*HA)
         PDUMA  =UA
         WA     =CB88*0.5D+00*HA/SQRT(GRDAA)
         PDUMXA =TWO*WA*GRADXA
         PDUMYA =TWO*WA*GRADYA
         PDUMZA =TWO*WA*GRADZA
C     ----- BETA CONTRIBTUTION ----
         RHO13B=ROB**(ONE/3.0D+00)
         RHOEB =RHO13B**4
         EB    =FACE
         IF(ROB.LT.1.0D-15) THEN
            GB    =ZERO
            HB    =ZERO
            PDUMB =ZERO
            PDUMXB=ZERO
            PDUMYB=ZERO
            PDUMZB=ZERO
            GOTO 9999
         ENDIF
         XB    =SQRT(GRDBB)/RHOEB
         QB    =SQRT(ONE+XB**2)
         SB    =LOG(XB+QB)
         TB    =ONE+6.0D+00*BECKEB*XB*SB
         GB    =-BECKEB*XB*XB/TB
         HB    =(   -TWO* BECKEB*XB
     >        +6.0D+00*(BECKEB*XB  )**2*(-SB+XB/QB))/TB**2
         RB    =FACP*RHO13B
         UB    =CB88*RB*(EB+GB-XB*HB)
         PDUMB  =UB
         WB    =CB88*0.5D+00*HB/SQRT(GRDBB)
         PDUMXB =TWO*WB*GRADXB
         PDUMYB =TWO*WB*GRADYB
         PDUMZB =TWO*WB*GRADZB
 9999    CONTINUE
         XCLA  = EA*RHOEA
         XCLB  = EB*RHOEB
         XCNLA = GA*RHOEA
         XCNLB = GB*RHOEB
         XALFA = CB88*XCLA* FTOTWT
         XALFB = CB88*XCLB* FTOTWT
         B88XA = CB88*XCNLA*FTOTWT
         B88XB = CB88*XCNLB*FTOTWT
C     ----- FOR OP CORRELATION -----
         SIGKA= -(FACE+GA)*TWO
         SIGKB= -(FACE+GB)*TWO
         DKDXA= -HA*TWO
         DKDXB= -HB*TWO
C     -------------------------------
         IF (LCFLAG) THEN
            CALL LRCFUN(UROHF,ROA,ROB,GRDAA,GRDBB,
     >           GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >           FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >           XALFA,XALFB,B88XA,B88XB,PDUMA,PDUMXA,PDUMYA,
     >           PDUMZA,PDUMB,PDUMXB,PDUMYB,PDUMZB,FACE)
         END IF
         XALF   = XALF + XALFA + XALFB
         B88X   = B88X + B88XA + B88XB
         DUMA   = DUMA + PDUMA
         DUMXA  = DUMXA+ PDUMXA
         DUMYA  = DUMYA+ PDUMYA
         DUMZA  = DUMZA+ PDUMZA
         DUMB   = DUMB + PDUMB
         DUMXB  = DUMXB+ PDUMXB
         DUMYB  = DUMYB+ PDUMYB
         DUMZB  = DUMZB+ PDUMZB
       ENDIF
      RETURN
      END
C*MODULE DFTXCA  *DECK LRCFUN
      SUBROUTINE LRCFUN(UROHF,ROA,ROB,GRDAA,GRDBB,
     >           GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >           FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >           XALFA,XALFB,B88XA,B88XB,DUMA,DUMXA,DUMYA,
     >           DUMZA,DUMB,DUMXB,DUMYB,DUMZB,FACE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION KS(2),MS(2),JUDGE
      LOGICAL UROHF,LCFLAG,LRINT
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      INTEGER LRFILE,NSPIN,S,I,P,J,A,IRS,IGSS,ITS,ILRC,NADD,MPROD,
     *        MAXPRD
      PARAMETER (NADD=1,MPROD=3,MAXPRD=3,NVAR=3)
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00,TWO=2.0D+00,THREE=3.0D+00,
     *           FOUR=4.0D+00,SIX=6.0D+00,EIGHT=8.0D+00,HALF=0.5D+00)
      PARAMETER (IRS=1,IGSS=2,ITS=3)
      PARAMETER (TAYTHR=1.0D-02,CUTLRC=1.0D-15,CUTHF=1.0D-15,
     >            CUTOFF=1.0D-15)
      DIMENSION FIPS(NADD,MAXPRD,2),DFIPS(NADD,MAXPRD,NVAR,2),
     *          PRDFIP(NADD,2),GS(2),ILRC(2),AS(2),EXPS1(2),
     *          DEXS(3,2),RS(2),DKSDXS(2),GFIPS(2),DRS(NVAR,2),
     *          XS(2),DXS(NVAR,2),DAS(NVAR,2),DMSDAS(2),DLSDAS(2),
     *          PRXALF(2),PRB88X(2)
C
C     INITIALIZE
C     NSPIN RHF=1 UHF=2
      IF (UROHF) THEN
         NSPIN = 2
      ELSE
         NSPIN = 1
      END IF
C-------------------------------
C     GET FIPS(FOR ENERGY)
C-------------------------------
      FACP     = FOUR/THREE
      ONETHR   = ONE/THREE
      PI12     = SQRT(PI)
      RS(1)    = ROA
      RS(2)    = ROB
      KS(1)    = SIGKA
      KS(2)    = SIGKB
C     CLFIPS AND GFIPS ARE DIVINDING PARAMETERS' FOR XALF AND B88X
      CLFIPS    = -TWO*FACE
      GS(1)    = GRDAA
      GS(2)    = GRDBB
      DKSDXS(1)= DKDXA
      DKSDXS(2)= DKDXB
      ILRC(1)  = 1
      ILRC(2)  = 1
C
      DO 10 S=1,NSPIN
         FIPS(1,1,S)=RS(S)**FACP
         FIPS(1,2,S)=KS(S)
         GFIPS(S)=KS(S)-CLFIPS
         AS(S)=(RS(S)**(-ONETHR))*(FIPS(1,2,S)**HALF)*EMU/(SIX*PI12)
         ERFVAL=DFTERF(HALF/AS(S))
         X4A21 =-ONE/(FOUR*(AS(S)**2))
         IF(ABS(X4A21).LT.TAYTHR) THEN
            EXPS1(S)=TAYEXP(X4A21)
         ELSE
            EXPS1(S)=EXP(X4A21)-ONE
         END IF
         FIRST = PI12*ERFVAL
         SE1   = TWO*EXPS1(S)
         SE2   = -FOUR*(AS(S)**2)*EXPS1(S)-ONE
         SECOND= AS(S)*(SE1+SE2)
         MS(S) = FIRST + SECOND
         FIPS(1,3,S)=-HALF+FACP*AS(S)*MS(S)
C
         JUDGE = ABS(AS(S) * MS(S) * (EIGHT/THREE))
C     DON'T APPLY LRC:DFT EXCHANGE FUNCTIONAL
         IF(JUDGE .LT. CUTLRC) THEN
            ILRC(S) = 0
C     FULL HARTREE-FOCK
         ELSE IF((ONE - JUDGE) .LT. CUTHF) THEN
            ILRC(S) = -1
         END IF
 10   CONTINUE
C
C     GET_PRODFIPS
      DO 20 S=1,NSPIN
C     FOR DIVIDE XALF AND B88X
      PRXALF(S) = ZERO
      PRB88X(S) = ZERO
         IF(ILRC(S).EQ.1) THEN
            DO 30 I=1,NADD
               PRDFIP(I,S)=ONE
               DO 40 P=1,MPROD,2
                  PRDFIP(I,S) = PRDFIP(I,S)*FIPS(I,P,S)
 40            CONTINUE
               PRXALF(S) = CLFIPS    * PRDFIP(I,S)
               PRB88X(S) = GFIPS(S) * PRDFIP(I,S)
C     FOR DFIPS
               PRDFIP(I,S)=PRDFIP(I,S)*FIPS(I,2,S)
 30         CONTINUE
         END IF
 20   CONTINUE
C     TWO/S:IF RHF EXC MUST PRXALF *2
      IF(ILRC(1).NE.0) THEN
         XALFA = CB88*(TWO/NSPIN)*PRXALF(1) * FTOTWT
         B88XA = CB88*(TWO/NSPIN)*PRB88X(1) * FTOTWT
      END IF
      IF(ILRC(2).NE.0) THEN
         XALFB = CB88*PRXALF(2)*FTOTWT
         B88XB = CB88*PRB88X(2)*FTOTWT
      END IF
C-----------------------------------
C     GET DFIPS(FOR 1ST DERIVATIVE)
C-----------------------------------
      DO 70 S=1,NSPIN
         DRS(IRS,S)  = ONE
         DRS(IGSS,S) = ZERO
         DRS(ITS,S)  = ZERO
         DO 80 A=1,NVAR
            DFIPS(1,1,A,S)=FACP*FIPS(1,1,S)/RS(S)*DRS(A,S)
 80      CONTINUE
         XS(S)=(GS(S)**HALF)*(RS(S)**(-FACP))
         DXS(IRS,S)  = -FACP/RS(S)*XS(S)
         IF(NEXFG.EQ.1) THEN
            DXS(IGSS,S) = ZERO
         ELSE
            DXS(IGSS,S) = HALF / GS(S) * XS(S)
         END IF
         DXS(ITS,S)  = ZERO
         DO 90 A=1,NVAR
            DFIPS(1,2,A,S) = DKSDXS(S) * DXS(A,S)
 90      CONTINUE
         DO 100 A=1,NVAR
            DAS(A,S) = AS(S)*(-ONETHR/RS(S)*DRS(A,S)+HALF
     >           / FIPS(1,2,S)*DFIPS(1,2,A,S))
 100     CONTINUE
         DMSDAS(S) = -THREE-12.0D+00*(AS(S)**2)*EXPS1(S)
         DLSDAS(S) = FACP*(MS(S)+AS(S)*DMSDAS(S))
         DO 110 A=1,NVAR
            DFIPS(1,3,A,S) = DLSDAS(S)*DAS(A,S)
 110     CONTINUE
 70   CONTINUE
C-----------------------------------
C     GET FIRST DERIVATIVE(DEXS)
C-----------------------------------
      DO 120 S=1,NSPIN
         DO 130 A=1,NVAR
            DEXS(A,S) = ZERO
            IF(ILRC(S).EQ.1) THEN
               DO 140 I=1,NADD
                  DO 150 J=1,MPROD
                     IF(ABS(DFIPS(I,J,A,S)).GT.CUTOFF) THEN
C     CALCULATE PRD_TERM=FIPS(1)*FIPS(2)*FIPS(3)/FIPS(S)
                        IF(ABS(FIPS(I,J,S)).GT.CUTOFF) THEN
                           PRDTER = PRDFIP(I,S)/FIPS(I,J,S)
                        ELSE
                           PRDTER = ONE
                           DO 160 P=1,MPROD
                              IF (P.NE.J) THEN
                                 PRDTER=PRDTER*FIPS(I,P,S)
                              END IF
 160                       CONTINUE
                        END IF
C     PRD_TERM TERMINATION
                        DEXS(A,S) = DEXS(A,S)+DFIPS(I,J,A,S)*PRDTER
                     END IF
 150              CONTINUE
 140           CONTINUE
            END IF
 130     CONTINUE
 120  CONTINUE
C------------------------------------
C     SET FIRST DERIVATIVE (DUM**)
C------------------------------------
      IF(ILRC(1).NE.0) THEN
         DUMA = CB88*DEXS(1,1)
         DUMXA= CB88*TWO*DEXS(2,1)*GRADXA
         DUMYA= CB88*TWO*DEXS(2,1)*GRADYA
         DUMZA= CB88*TWO*DEXS(2,1)*GRADZA
      END IF
      IF(ILRC(2).NE.0) THEN
         IF(UROHF) THEN
            DUMB = CB88*DEXS(1,2)
            DUMXB= CB88*TWO*DEXS(2,2)*GRADXB
            DUMYB= CB88*TWO*DEXS(2,2)*GRADYB
            DUMZB= CB88*TWO*DEXS(2,2)*GRADZB
         END IF
      END IF
C     IF(ILRC.EQ.0) XALF,B88X,DUM** IS ORIGINAL DFT EXCHANGE FUNCTIONAL
      RETURN
      END
C
C*MODULE DFTXCA  *DECK TAYEXP
      DOUBLE PRECISION FUNCTION TAYEXP(X)
      DOUBLE PRECISION X
      INTEGER NTERM,I
      DOUBLE PRECISION SUMVAL,IVAL,COEF
      PARAMETER (NTERM=16)
C
      SUMVAL = X
      IVAL = X
      DO 10 I = 2,NTERM
         COEF =  I
         IVAL = IVAL * (X / COEF)
         SUMVAL = SUMVAL + IVAL
 10         CONTINUE
      TAYEXP = SUMVAL
      RETURN
      END
C
C*MODULE DFTXCA  *DECK DFTERF
      DOUBLE PRECISION FUNCTION DFTERF(X)
C
      DOUBLE PRECISION X
      DOUBLE PRECISION A,X2,Y,TMP,SER,COF,STP,GLN,EPS,FPMIN,
     *       GAMSER,AP,SUM,DEL,B,C,D,H,AN
      INTEGER J,ITMAX
      DIMENSION COF(6)
C
      DATA COF,STP/76.18009172947146D+00,-86.50532032941677D+00,
     *       24.01409824083091D+00,-1.231739572450155D+00,
     *       .1208650973866179D-02,-.5395239384953D-05,
     *       2.5066282746310005D+00/
C
C          THIS CALCULATES THE STANDARD ERF VALUE.
C          TO AVOID PROBLEMS WITH MATH LIBRARIES, GIVE IT A UNIQUE NAME.
C          CAN IT BE MORE RAPIDLY EVALUATED WITH GRAHAM FLETCHER'S CODE?
C
        ITMAX=100
        EPS=3.0D-07
        FPMIN=1.0D-30
        A=0.5D+00
        X2=X*X
C     GLN
        Y=A
        TMP=A+5.5D+00
        TMP=(A+0.5D+00)*LOG(TMP)-TMP
        SER=1.00000000190015D+00
        DO 30 J=1,6
           Y=Y+1.0D+00
           SER=SER+COF(J)/Y
 30             CONTINUE
        GLN=TMP+LOG(STP*SER/A)
C     GAMSER
        IF(X2.LT.A+1.0D+00) THEN
           IF(X2.LE.0.0D+00) THEN
              WRITE(6,9100)
              GAMSER = 0.0D+00
           ELSE
              AP = A
              SUM = 1.0D+00/A
              DEL = SUM
              DO 10 J=1,ITMAX
                 AP=AP+1.0D+00
                 DEL=DEL*X2/AP
                 SUM=SUM+DEL
                 IF(ABS(DEL).LT.ABS(SUM)*EPS) GOTO 100
 10                         CONTINUE
              WRITE(6,9200)
C
 100                    GAMSER=SUM*EXP(-X2+A*LOG(X2)-GLN)
           END IF
        ELSE
           B=X2+1.0D+00-A
           C=1.0D+00/FPMIN
           D=1.0D+00/B
           H=D
           DO 20 J=1,ITMAX
              AN=-J*(J-A)
              B=B+2.0D+00
              D=AN*D+B
              IF(ABS(D).LT.FPMIN) D=FPMIN
              C=B+AN/C
              IF(ABS(C).LT.FPMIN) C=FPMIN
              D=1.0D+00/D
              DEL=D*C
              H=H*DEL
              IF (ABS(DEL-1.0D+00).LT.EPS) GOTO 110
 20                   CONTINUE
           WRITE(6,9300)
C
 110              GAMSER=EXP(-X2+A*LOG(X2)-GLN)*H
           GAMSER = 1.0D+00 - GAMSER
        ENDIF
C
        IF(X.LT.0.0D+00) THEN
           DFTERF = -GAMSER
        ELSE
           DFTERF =  GAMSER
        END IF
C     ERF MUST LESS THAN 1.0
        IF (DFTERF.GT.1.0D+00) THEN
           WRITE (6,9400) DFTERF
        END IF
C
 9100      FORMAT('X < 0 IN GSER')
 9200      FORMAT('A TOO LARGE, ITMAX TOO SMALL IN GSER')
 9300      FORMAT('A TOO LARGE, ITMAX TOO SMALL IN GCF')
 9400      FORMAT('DFTERF VALUE > 1.0  DFTERF= ',E20.10)
        RETURN
        END
C*MODULE DFTXCA  *DECK CALCEXC
      SUBROUTINE CALCEXC(ROA,ROB,FTOTWT,GRDAA,GRDBB,GRDAB,
     >                  GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                  XALPHA,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >                  VXCB1,DUMBX,DUMBY,DUMBZ,ECF,
     >                  TAUXA,TAUYA,TAUZA,DMGGA,TAUXB,TAUYB,TAUZB,DMGGB)
C
C***********************************************************************
C       DRIVER FOR THE  DFT EXCHANGE
C      CORRELATION ENERGY CALCULATION
C***********************************************************************
C
C     NEXFG CONTROLS EXCHANGE FUNCTIONAL
C    -1     ---  BECKE 3 PARAMETER
C     0     ---  HARTREE-FOCK
C     1     ---  SLATER
C     2     ---  BECKE 88
C     3     ---  PERDEW-WANG 91
C     4     ---  PERDEW-BURKE-ERNZERHOF (PBE)
C     5     ---  GILL 96
C     6     ---  PARAMETER-FREE
C     7     ---  HANDY-COHEN (OPTX) NOT RELATED TO (OP) CORRELATION
C     11    ---  M05
C     12    ---  VS98
C     13    ---  M06
C     8,14  ...  PLEASE ADD NEW FUNCTIONAL
C
C     NCORFG CONTROLS CORRELATION FUNCTIONAL
C    -3     ---  PERDEW-WANG 91 USING NWCHEM CONSTANTS FOR PW91LDA
C    -1     ---  B3LYP
C     0     ---  NONE (EXCHANGE ONLY)
C     1     ---  VWN FORMULA 5
C     2     ---  LEE-YANG-PARR
C     3     ---  PERDEW-WANG 91
C     4     ---  PERDEW-BURKE-ERNZERHOF (PBE)
C     5     ---  VWN FORMULA 1 (RPA)
C     6     ---  VWN FORMULA 5 FOR BACKWARDS COMPATIBILITY
C     7     ---  ONE-PARAMETER-PROGRESSIVE (OP)
C     8     ---  PERDEW-WANG 91 LDA
C     9     ---  PERDEW 86
C     10    ---  PERDEW-ZUNGER 81
C     11    ---  M05
C     12    ---  VS98
C     13    ---  M06
C     14    ...  PLEASE ADD NEW FUNCTIONAL
C
C   * OTHER VARIABLES....
C
C     ROA,ROB                    <---     DENSITY
C     GRADXA,GRADXB,ETC          <---     GRAD(RHO)
C     GRDAA,GRDBB,GRDAB          <---     GRAD(RHO) * GRAD(RHO)
C     VXCA1,VXCB1                <---     DEXC/DRHO
C     DUMAX,DUMAY,DUMAZ,ETC      <---     DEXC/DG
C     XALPHA                     <---     EX FROM LDA TERM
C     XGRD                       <---     EX FROM GRADIENT TERM
C     ECF                        <---     CORRELATION ENERGY
C     TAU_A, TAU_B               <---     X,Y,Z COMPONENTS OF KINETIC
C                                         ENERGY DENSITY
C     DMGGA, DMGGB               <---     DEXC/DTAU
C***********************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /CM05  / CM05XF, IM05
      COMMON /CM06  / IM06
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      LOGICAL UROHF, RHFGVB
C
      DATA RHF,UHF,ROHF,GVB/8HRHF     ,8HUHF     ,8HROHF     ,
     *     8HGVB     /
C***********************************************************************
      UROHF = SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
      RHFGVB = SCFTYP.EQ.RHF  .OR.  SCFTYP.EQ.GVB
C
C***********************************************************************
C           EXCHANGE FUNCTIONAL CALCULATION
C***********************************************************************
      DKDXA=0.0D+00
      DKDXB=0.0D+00
      SIGKA=0.0D+00
      SIGKB=0.0D+00
      DUMMY=0.0D+00
C
C     -----HARTREE-FOCK EXCHANGE-----
C
      IF(NEXFG.EQ.0) THEN
         IF (NCORFG.EQ.7) THEN
C     ----- HOP METHOD -----
            CALL B88XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >           GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >           FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >           DUMMY,DUMMY,DUMMY,DUMMY,DUMMY,DUMMY,
     >           DUMMY,DUMMY,DUMMY,DUMMY)
C
         ENDIF
C
C     ----- SLATER EXCHANGE -----
C
      ELSE IF(NEXFG.EQ.1) THEN
         CALL LSDSLT(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA1,VXCB1,XALPHA,
     *               SIGKA,SIGKB)
C
C     ----- BECKE 88 EXCHANGE -----
C
      ELSE IF(NEXFG.EQ.2) THEN
C
         CALL B88XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >        XALPHA,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >        VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- PW91 EXCHANGE -----
C
      ELSE IF(NEXFG.EQ.3) THEN
      CALL PW91XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >        VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- PBE EXCHANGE -----
C
C     CAUTION: THIS FUNCTIONAL MAY NOT BE CORRECTLY PROGRAMMED.
C     MIKE SCHMIDT FOUND IN JUNE 2001 THAT ETHANOL 6-31G(D) HAS
C     AN ENERGY OF -153.9850674001 WITH GAMESS,
C              BUT -153.9851024042 WITH NWCHEM.  SINCE THE OTHER
C     FUNCTIONALS AGREE TO 1 OR 2 IN THE 7TH DIGIT, THIS SEEMS
C     TO BE A FAIRLY LARGE DISCREPANCY.  NWCHEM RUN WITH "XPBE".
C
C         PBE EXCHANGE/CORRELATION FUNCTIONAL:
C           J.P.PERDEW, K.BURKE, M.ERNZERHOF
C              PHYS.REV.LETT.  77, 3865-8(1996); ERR. 78,1396(1997)
C           M.ERNZERHOF, G.E.SCUSERIA
C              J.CHEM.PHYS. 110, 5029-5036(1999)
C         NOTE THAT ONLY THE EXCHANGE FUNCTIONAL IS IN GAMESS.
C
C     ISSUE ADDRESSED BY SAROM SOK ON 03/02/2008
C     THE VALUE FOR MU NOW HAS INCREASE PRECISION TO MATCH OTHER CODES
C
      ELSE IF(NEXFG.EQ.4) THEN
         CALL PBEXFOP(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >        XALPHA,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >        VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- GILL 96 EXCHANGE -----
C
      ELSE IF(NEXFG.EQ.5) THEN
         CALL G96XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >        XALPHA,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >        VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- PARAMETER FREE EXCHANGE -----
C
      ELSE IF(NEXFG.EQ.6) THEN
         CALL PFREE(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >        XALPHA,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >        VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- OPTX EXCHANGE -----
C
      ELSE IF(NEXFG.EQ.7) THEN
         CALL LSDSLT(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA1,VXCB1,XALPHA,
     >               SIGKA,SIGKB)
         CALL OPTXXF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >               GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >               FTOTWT,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >               VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- M05 EXCHANGE -----
C
      ELSE IF(NEXFG.EQ.11) THEN
         CALL M05XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,XGRD,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,DMGGA,
     >        VXCB1,DUMBX,DUMBY,DUMBZ,DMGGB,
     >        TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IM05)
C
C     ----- VS98 EXCHANGE -----
C
      ELSE IF(NEXFG.EQ.12) THEN
         CALL VS98XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,XGRD,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,DMGGA,
     >        VXCB1,DUMBX,DUMBY,DUMBZ,DMGGB,
     >        TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IM06)
C
C     ----- M06 EXCHANGE -----
C
      ELSE IF(NEXFG.EQ.13) THEN
         CALL M06XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,XGRD,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,DMGGA,
     >        VXCB1,DUMBX,DUMBY,DUMBZ,DMGGB,
     >        TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IM06)
C
C     ----- B3 EXCHANGE -----
C
      ELSE IF(NEXFG.EQ.-1) THEN
         CALL LSDSLT(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA1,VXCB1,XALPHA,
     >               SIGKA,SIGKB)
         CALL B88XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >        XALPHA,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >        VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- X3LYP EXCHANGE -----
C
       ELSE IF(NEXFG.EQ.-2) THEN
         CALL LSDSLT(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA1,VXCB1,XALPHA,
     >        SIGKA,SIGKB)
         CALL B88XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >        XALPHA,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >        VXCB1,DUMBX,DUMBY,DUMBZ)
         CALL PW91XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >        VXCB1,DUMBX,DUMBY,DUMBZ)
      ELSE
         WRITE(IW,*) 'CANNOT USE THIS EXCHANGE FUNCTIONAL'
         CALL ABRT
         STOP
      ENDIF
C***********************************************************************
C           CORRELATION FUNCTIONAL CALCULATION
C***********************************************************************
C
C     ----- NO CORRELATION FUNCTIONAL-----
C
      IF (NCORFG.EQ.0) THEN
C
C     ----- VWN CORRELATION FUNCTIONAL -----
C
      ELSE IF(NCORFG.EQ.1) THEN
         CALL VWN5CF(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA1,VXCB1,ECF)
C
C     ----- LEE-YANG-PARR CORRELATION FUNCTIONAL -----
C
      ELSE IF(NCORFG.EQ.2) THEN
         CALL LYPCF(RHFGVB,UROHF,NB,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,VXCB1,DUMBX,DUMBY,DUMBZ)
C
C
C     ----- PERDEW-WANG 91 CORRELATION FUNCTIONAL-----
C
      ELSE IF(NCORFG.EQ.3) THEN
         CALL PW91CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- PBE CORRELATION FUNCTIONAL-----
C
      ELSE IF(NCORFG.EQ.4) THEN
         CALL PBECF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- VWN FORMULA 1 CORRELATION FUNCTIONAL -----
C
      ELSE IF(NCORFG.EQ.5) THEN
         CALL VWN1CF(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA1,VXCB1,ECF)
C
C     ----- ONE-PARAMETER PROGRESSIVE CORRELATION FUNCTIONAL -----
C
      ELSE IF(NCORFG.EQ.7) THEN
         CALL OPCOR(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        SIGKA,SIGKB,DKDXA,DKDXB,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- PERDEW-WANG 91 LDA CORRELATION FUNCTIONAL -----
C
      ELSE IF(NCORFG.EQ.8) THEN
         CALL PW91LDACF(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA1,VXCB1,ECF)
C
C     ----- PERDEW 86 CORRELATION FUNCTIONAL -----
C
      ELSE IF(NCORFG.EQ.9) THEN
         CALL P86CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- PERDEW ZUNGER 81 CORRELATION FUNCTIONAL -----
C
      ELSE IF(NCORFG.EQ.10) THEN
         CALL PZ81CF(RHFGVB,UROHF,ROA,ROB,FTOTWT,ECF,VXCA1,VXCB1)
C
C     ----- M05 CORRELATION FUNCTIONAL -----
C
      ELSE IF(NCORFG.EQ.11) THEN
         CALL M05CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,DMGGA,
     >        VXCB1,DUMBX,DUMBY,DUMBZ,DMGGB,
     >        TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IM05)
C
C     ----- VS98 CORRELATION FUNCTIONAL -----
C
      ELSE IF(NCORFG.EQ.12) THEN
         CALL VS98CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,DMGGA,
     >        VXCB1,DUMBX,DUMBY,DUMBZ,DMGGB,
     >        TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IM06)
C
C     ----- M06 CORRELATION FUNCTIONAL -----
C
      ELSE IF(NCORFG.EQ.13) THEN
         CALL M06CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,DMGGA,
     >        VXCB1,DUMBX,DUMBY,DUMBZ,DMGGB,
     >        TAUXA,TAUYA,TAUZA,TAUXB,TAUYB,TAUZB,IM06)
C
C     ----- B3LYP CORRELATION PART -----
C
      ELSE IF(NCORFG.EQ.-1) THEN
         CALL VWN5CF(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA1,VXCB1,ECF)
         CALL LYPCF(RHFGVB,UROHF,NB,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- GAUSSIAN 94 B3LYP CORRELATION PART -----
C     -----             X3LYP CORRELATION PART -----
C
      ELSE IF(NCORFG.EQ.-2) THEN
         CALL VWN1CF(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA1,VXCB1,ECF)
         CALL LYPCF(RHFGVB,UROHF,NB,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,VXCB1,DUMBX,DUMBY,DUMBZ)
C
C     ----- PERDEW-WANG 91 CORRELATION FUNCTIONAL -----
C     ----- UTILIIZING NWCHEM CONSTANTS FOR PW91LDA ---
C
      ELSE IF(NCORFG.EQ.-3) THEN
         CALL NWPW91CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,VXCA1,DUMAX,DUMAY,DUMAZ,
     >        VXCB1,DUMBX,DUMBY,DUMBZ)
C
      ELSE
         WRITE(IW,*) 'CANNOT USE THIS CORRELATIONAL FUNCTIONAL'
         CALL ABRT
         STOP
      ENDIF
      RETURN
      END
C*MODULE DFTXCA  *DECK G96XF
      SUBROUTINE G96XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                 GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                 FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >                 XALF,B88X,DUMA,DUMXA,DUMYA,DUMZA,
     >                 DUMB,DUMXB,DUMYB,DUMZB)
C***********************************************************************
C     GILL 96 EXCHANGE FUNCTIONAL
C
C     P. M. W. GILL
C     MOL. PHYS. 89, 433 (1996)
C
C     -----
C     INPUT
C     -----
C
C     ROA,B   .... SPIN DENSITIES
C     GRAD**  .... DENSITY GRADIENTS
C     TOTWT   .... VOLUME OF THIS GRID
C
C     ------
C     OUTPUT
C     ------
C
C     XALF    ....  EXCHANGE ENERGY OF LDA PART
C     B88X    ....  EXCHANGE ENERGY OF GGA PART
C     DUM***  ....  FIRST DERIVATIVES OF THE EXCHANGE FUNCTINALS
C     SIGKA,B ....  EXC=-1/2*INT(SIGK*RHO^(4/3))
C     DKDXA,B ....  DSIGK/DX
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00)
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      LOGICAL LCFLAG,LRINT
C***********************************************************************
       X1 = (3.0D+00/PI)**(ONE/3.0D+00)
       TWO13 = TWO**(ONE/3.0D+00)
       FACE  =-0.75D+00*X1*TWO13
       FACP  = 4.0D+00/3.0D+00
C***********************************************************************
C      CLOSED SHELL CALCULATION
C***********************************************************************
       IF (RHFGVB) THEN
         RHO=ROA
         GRD=GRDAA
         RHO13=RHO**(ONE/3.0D+00)
         RHOE =RHO13*RHO
         X    =SQRT(GRD)/RHOE
         G    =-X**(THREE/TWO)/137.0D+00
         H    =-(1.5D+00/137.0D+00)*SQRT(X)
         R    =FACP*RHO13
         E    =FACE
         XCL  =E*RHOE
         XCNL =G*RHOE
         XALFA =CSLT*XCL *TWO*FTOTWT
         B88XA =CB88*XCNL*TWO*FTOTWT
         U     =CSLT*R*(E+G-X*H)
         PDUMA =U
         W     =CB88*0.5D+00*H/SQRT(GRD)
         PDUMXA=TWO*W*GRADXA
         PDUMYA=TWO*W*GRADYA
         PDUMZA=TWO*W*GRADZA
C     ----- FOR OP CORRELATION -----
         SIGKA= -(FACE+G)*TWO
         DKDXA= -H*TWO
C     ------------------------------
         IF (LCFLAG) THEN
            CALL LRCFUN(UROHF,ROA,ROB,GRDAA,GRDBB,
     >           GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >           FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >           XALFA,XALFB,B88XA,B88XB,PDUMA,PDUMXA,PDUMYA,
     >           PDUMZA,PDUMB,PDUMXB,PDUMYB,PDUMZB,FACE)
         END IF
         XALF   = XALF + XALFA
         B88X   = B88X + B88XA
         DUMA   = DUMA + PDUMA
         DUMXA  = DUMXA+ PDUMXA
         DUMYA  = DUMYA+ PDUMYA
         DUMZA  = DUMZA+ PDUMZA
C***********************************************************************
C      OPEN SHELL CALCULATION
C***********************************************************************
       ELSE IF (UROHF) THEN
         RHO13A=ROA**(ONE/3.0D+00)
         RHOEA =RHO13A*ROA
         EA    =FACE
         IF(ROA.LT.1.0D-10) THEN
           GA  = ZERO
           PDUMA=ZERO
           PDUMXA=ZERO
           PDUMYA=ZERO
           PDUMZA=ZERO
           GOTO 9998
         ENDIF
         XA     =SQRT(GRDAA)/RHOEA
         GA     =-XA**(THREE/TWO)/137.0D+00
         HA     =-(1.5D+00/137.0D+00)*SQRT(XA)
         RA     =FACP*RHO13A
         UA     =CSLT*RA*(EA+GA-XA*HA)
         PDUMA  =UA
         WA     =CB88*0.5D+00*HA/SQRT(GRDAA)
         PDUMXA =TWO*WA*GRADXA
         PDUMYA =TWO*WA*GRADYA
         PDUMZA =TWO*WA*GRADZA
 9998    CONTINUE
         RHO13B=ROB**(ONE/3.0D+00)
         RHOEB =RHO13B*ROB
         EB    =FACE
         IF(ROB.LT.1.0D-10) THEN
           GB  = ZERO
           PDUMB=ZERO
           PDUMXB=ZERO
           PDUMYB=ZERO
           PDUMZB=ZERO
           GOTO 9999
         ENDIF
         XB    =SQRT(GRDBB)/RHOEB
         GB    =-XB**(THREE/TWO)/137.0D+00
         HB    =-(1.5D+00/137.0D+00)*SQRT(XB)
         RB    =FACP*RHO13B
         UB    =CSLT*RB*(EB+GB-XB*HB)
         PDUMB =UB
         WB    =CB88*0.5D+00*HB/SQRT(GRDBB)
         PDUMXB =TWO*WB*GRADXB
         PDUMYB =TWO*WB*GRADYB
         PDUMZB =TWO*WB*GRADZB
 9999    CONTINUE
         XCLA   = EA*RHOEA
         XCLB   = EB*RHOEB
         XCNLA  = GA*RHOEA
         XCNLB  = GB*RHOEB
         XALFA  = CSLT*XCLA *FTOTWT
         XALFB  = CSLT*XCLB *FTOTWT
         B88XA  = CB88*XCNLA*FTOTWT
         B88XB  = CB88*XCNLB*FTOTWT
C      ----- FOR OP CORRELATION -----
         SIGKA= -(FACE+GA)*TWO
         SIGKB= -(FACE+GB)*TWO
         DKDXA= -HA*TWO
         DKDXB= -HB*TWO
C     -------------------------------
         IF (LCFLAG) THEN
            CALL LRCFUN(UROHF,ROA,ROB,GRDAA,GRDBB,
     >           GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >           FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >           XALFA,XALFB,B88XA,B88XB,PDUMA,PDUMXA,PDUMYA,
     >           PDUMZA,PDUMB,PDUMXB,PDUMYB,PDUMZB,FACE)
         END IF
         XALF   = XALF + XALFA + XALFB
         B88X   = B88X + B88XA + B88XB
         DUMA   = DUMA + PDUMA
         DUMXA  = DUMXA+ PDUMXA
         DUMYA  = DUMYA+ PDUMYA
         DUMZA  = DUMZA+ PDUMZA
         DUMB   = DUMB + PDUMB
         DUMXB  = DUMXB+ PDUMXB
         DUMYB  = DUMYB+ PDUMYB
         DUMZB  = DUMZB+ PDUMZB
      ENDIF
      RETURN
      END
C*MODULE DFTXCA  *DECK INPGDFT
      SUBROUTINE INPGDFT(DFTYPE,PFTYP)
C***********************************************************************
C     NEXFG CONTROLS EXCHANGE FUNCTIONAL
C    -1     ---  BECKE 3 PARAMETER
C     0     ---  HARTREE-FOCK
C     1     ---  SLATER
C     2     ---  BECKE 88
C     3     ---  PERDEW-WANG 91
C     4     ---  PERDEW-BURKE-ERNZERHOF (PBE)
C     5     ---  GILL 96
C     6     ---  PARAMETER-FREE
C     7     ---  HANDY-COHEN (OPTX) NOT RELATED TO (OP) CORRELATION
C     11    ---  M05, M05-2X (CONTROLLED BY IM05)
C     12    ---  VS98
C     13    ---  M06, M06-L, M06-HF, M06-2X (CONTROLLED BY IM06)
C     8,14  ...  PLEASE ADD NEW FUNCTIONAL
C
C     NCORFG CONTROLS CORRELATION FUNCTIONAL
C    -3     ---  PERDEW-WANG 91 USING NWCHEM CONSTANTS FOR PW91LDA
C    -1     ---  B3LYP
C     0     ---  NONE (EXCHANGE ONLY)
C     1     ---  VWN FORMULA 5
C     2     ---  LEE-YANG-PARR
C     3     ---  PERDEW-WANG 91
C     4     ---  PERDEW-BURKE-ERNZERHOF (PBE)
C     5     ---  VWN FORMULA 1 (RPA)
C     6     ---  VWN FORMULA 5 FOR BACKWARDS COMPATIBILITY
C     7     ---  ONE-PARAMETER-PROGRESSIVE (OP)
C     8     ---  PERDEW-WANG 91 LDA
C     9     ---  PERDEW 86
C     10    ---  PERDEW-ZUNGER 81
C     11    ---  M05, M05-2X (CONTROLLED BY IM05)
C     12    ---  VS98
C     13    ---  M06, M06-L, M06-HF, M06-2X (CONTROLLED BY IM06)
C     14    ...  PLEASE ADD NEW FUNCTIONAL
C***********************************************************************
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C      /DFTHF/,/DFTCF/ ARE INFORMATION FOR GRID CODE
C      /DFTPAR/ IS INFORMATION FOR GRID-FREE CODE
C
      SAVE EXCHR,CORCHR
C
      LOGICAL LCFLAG,LRINT,NEEDTAU,PRTTAU,SG1
C
      PARAMETER (MXGRID=10)
C
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /DFLEB/  NLEB(MXGRID),NLEB0(MXGRID)
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /METGGA/ NEEDTAU,PRTTAU
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      COMMON /CM05  / CM05XF, IM05
      COMMON /CM06  / IM06
      COMMON /CPBE0 / CPBEX
      COMMON /CX3LYP/ CPW91
      COMMON /CPW91L/ CPW91LDA
      COMMON /EXFACT/ COPTX
C
C     NUMBER OF EXCHANGE AND CORRELATION FUNCTIONALS
C
      PARAMETER (NEX=10, NCOR=14)
      CHARACTER*8  EX(NEX),COR(NCOR)
      CHARACTER*10 EXNM(NEX),CORNM(NCOR)
      DIMENSION    NUMEX(NEX),NUMCOR(NCOR),DQOP(NEX)
C
      CHARACTER*10 EXCHR,CORCHR,XCHNG,CRRN
      LOGICAL GOPARR,DSKWRK,MASWRK
      DOUBLE PRECISION LOCAL,LSDA,LCBOP,LCBLYP
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
C     NOTE THAT 'PFREE' IS NOT A DOCUMENTED EXCHANGE FUNCTIONAL
C
      DATA EX/'S','B','PW91','PBE','G','P','O','M05','VS98','M06'/
      DATA NUMEX/1,2,3,4,5,6,7,11,12,13/
      DATA EXNM/'SLATER','BECKE88','PW91','PBE','GILL96','PFREE','OPTX',
     >          'M05','VS98','M06'/
      DATA DQOP/2.5654D+00,2.3670D+00,ZERO,2.3789D+00,2.4167D+00,
     *          2.3866D+00,ZERO,ZERO,ZERO,ZERO/
C
C     NOTE THAT VWN5 IS NOT A DOCUMENTED SYNONYM FOR VWN.
C     NOTE THAT PW91LDA IS NOT A DOCUMENTED CORRELATION FUNCTIONAL.
C     NOTE THAT NW91 IS NOT A DOCUMENTED CORRELATION FUNCTIONAL.
C
      DATA COR/'VWN5','LYP','PW91','PBE','VWN1','VWN','OP',
     +           'PW91LDA','P86','PZ81','NW91','M05','VS98','M06'/
      DATA NUMCOR/1,2,3,4,5,6,7,8,9,10,-3,11,12,13/
      DATA CORNM/'VWN5','LYP88','PW91','PBE','VWN1','VWN','OP',
     +           'PW91LDA','P86','PZ81','PW91-NW','M05','VS98','M06'/
C
      DATA LOCAL    /8HLOCAL   /,SLATER /8HSLATER  /,BECKE  /8HBECKE   /
      DATA GILL     /8HGILL    /,G96    /8HG96     /
      DATA HOP      /8HOP      /,HLYP   /8HLYP     /
      DATA LSDA     /8HLSDA    /,PFREE  /8HPFREE   /
      DATA B3LYP    /8HB3LYP   /,BHHLYP /8HBHHLYP  /,HALF   /8HHALF    /
      DATA LCBOP    /8HLCBOP   /,LCBLYP /8HLCBLYP  /
      DATA HVWN1    /8HVWN1    /,HVWN5  /8HVWN5    /,HVWN   /8HVWN     /
      DATA PW91     /8HPW91    /,PW91X  /8HPW91X   /,HPW91  /8HPW91C   /
      DATA PBE      /8HPBE     /,PBEX   /8HPBEX    /,HPBE   /8HPBEC    /
      DATA PBE0     /8HPBE0    /,PBE0M  /8HPBE0M   /
      DATA B3LYP1   /8HB3LYP1  /,X3LYP  /8HX3LYP   /,OPTX   /8HOPTX    /
      DATA HPW91LDA /8HPW91LDA /,P86    /8HP86     /,PZ81   /8HPZ81    /
      DATA HPW91NW  /8HNW91    /
C
C     M06 KEYWORDS
C
      DATA  TRUHLAR05   /8HM05     /
      DATA CTRUHLAR05   /8HM05C    /
      DATA XTRUHLAR05   /8HM05X    /
      DATA  TRUHLAR052X /8HM05-2X  /
      DATA CTRUHLAR052X /8HM05-2XC /
      DATA XTRUHLAR052X /8HM05-2XX /
      DATA  TRUHLAR06   /8HM06     /
      DATA CTRUHLAR06   /8HM06C    /
      DATA XTRUHLAR06   /8HM06X    /
      DATA  TRUHLAR06L  /8HM06-L   /
      DATA CTRUHLAR06L  /8HM06-LC  /
      DATA XTRUHLAR06L  /8HM06-LX  /
      DATA  TRUHLAR062X /8HM06-2X  /
      DATA CTRUHLAR062X /8HM06-2XC /
      DATA XTRUHLAR062X /8HM06-2XX /
      DATA  TRUHLAR06HF /8HM06-HF  /
      DATA CTRUHLAR06HF /8HM06-HFC /
      DATA XTRUHLAR06HF /8HM06-HFX /
      DATA  VS98        /8HVS98    /
      DATA CVS98        /8HVS98C   /
      DATA XVS98        /8HVS98X   /
C
      DATA CHECK/8HCHECK   /
C
      NDFTFG=1
      NEXFG =0
      NCORFG=0
      CSLT     = ONE
      CB88     = ONE
      CLYP     = ONE
      CVWN     = ONE
      CPW91    = ONE
      CPBEX    = ONE
      CPW91LDA = ONE
      CM05XF   = ONE
      COPTX    = 1.43169D+00
      DFTTYP(3)= 0.0D+00
C
C     ----- EXCHANGE ONLY METHOD -----
C
      IF(DFTYPE.EQ.LOCAL .OR. DFTYPE.EQ.LSDA .OR. DFTYPE.EQ.SLATER) THEN
         NEXFG     = 1
         EXCHR     = EXNM(1)
         CORCHR    = 'NONE'
      ELSEIF(DFTYPE.EQ.BECKE)  THEN
         NEXFG     = 2
         EXCHR     = EXNM(2)
         CORCHR    = 'NONE'
      ELSEIF(DFTYPE.EQ.PW91X) THEN
         NEXFG     = 3
         EXCHR     = EXNM(3)
         CORCHR    = 'NONE'
      ELSEIF(DFTYPE.EQ.PBEX)  THEN
         NEXFG     = 4
         EXCHR     = EXNM(4)
         CORCHR    = 'NONE'
      ELSEIF(DFTYPE.EQ.G96.OR.DFTYPE.EQ.GILL)  THEN
         NEXFG     = 5
         EXCHR     = EXNM(5)
         CORCHR    = 'NONE'
      ELSEIF(DFTYPE.EQ.PFREE)  THEN
         NEXFG     = 6
         EXCHR     = EXNM(6)
         CALL INPPFREE(PFTYP)
         CORCHR    = 'NONE'
      ELSEIF(DFTYPE.EQ.OPTX)  THEN
         NEXFG     = 7
         EXCHR     = EXNM(7)
         CORCHR    = 'NONE'
         CSLT      = 1.05151D+00
      ELSEIF(DFTYPE.EQ.XTRUHLAR05) THEN
         NEXFG     = 11
         EXCHR     = EXNM(8)
         CORCHR    = 'NONE'
         IM05      = 1
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.XTRUHLAR052X) THEN
         NEXFG     = 11
         EXCHR     = EXNM(8)
         CORCHR    = 'NONE'
         IM05      = 2
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.XVS98) THEN
         NEXFG     = 12
         EXCHR     = EXNM(9)
         CORCHR    = 'NONE'
         IM06      = 1
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.XTRUHLAR06) THEN
         NEXFG     = 13
         EXCHR     = EXNM(10)
         CORCHR    = 'NONE'
         IM06      = 3
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.XTRUHLAR06L) THEN
         NEXFG     = 13
         EXCHR     = 'M06-L'
         CORCHR    = 'NONE'
         IM06      = 1
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.XTRUHLAR062X) THEN
         NEXFG     = 13
         EXCHR     = 'M06-2X'
         CORCHR    = 'NONE'
         IM06      = 4
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.XTRUHLAR06HF) THEN
         NEXFG     = 13
         EXCHR     = 'M06-HF'
         CORCHR    = 'NONE'
         IM06      = 2
         NEEDTAU   = .TRUE.
C
C    ----- CORRELATION ONLY METHOD -----
C
      ELSEIF(DFTYPE.EQ.HVWN)   THEN
         NCORFG    = 1
         EXCHR     = 'HFX'
         CORCHR    = CORNM(1)
         DFTTYP(3) = 1.0D+00
      ELSEIF(DFTYPE.EQ.HVWN5)   THEN
         NCORFG    = 1
         EXCHR     = 'HFX'
         CORCHR    = CORNM(1)
         DFTTYP(3) = 1.0D+00
      ELSEIF(DFTYPE.EQ.HLYP)   THEN
         NCORFG    = 2
         EXCHR     = 'HFX'
         CORCHR    = CORNM(2)
         DFTTYP(3) = 1.0D+00
      ELSEIF(DFTYPE.EQ.HPW91)   THEN
         NCORFG    = 3
         EXCHR     = 'HFX'
         CORCHR    = CORNM(3)
         DFTTYP(3) = 1.0D+00
      ELSEIF(DFTYPE.EQ.HPW91NW)  THEN
         NCORFG    =-3
         EXCHR     = 'HFX'
         CORCHR    = 'PW91-NW'
         DFTTYP(3) = 1.0D+00
      ELSEIF(DFTYPE.EQ.HPBE)   THEN
         NCORFG    = 4
         EXCHR     = 'HFX'
         CORCHR    = CORNM(4)
         DFTTYP(3) = 1.0D+00
      ELSEIF(DFTYPE.EQ.HVWN1)   THEN
         NCORFG    = 5
         EXCHR     = 'HFX'
         CORCHR    = CORNM(5)
         DFTTYP(3) = 1.0D+00
      ELSEIF(DFTYPE.EQ.HOP)     THEN
         NCORFG    = 7
         EXCHR     = 'HFX'
         CORCHR    = CORNM(7)
         DFTTYP(3) = 1.0D+00
         QOP       = 2.3670D+00
      ELSEIF(DFTYPE.EQ.HPW91LDA) THEN
         NCORFG    = 8
         EXCHR     = 'HFX'
         CORCHR    = CORNM(8)
         DFTTYP(3) = 1.0D+00
      ELSEIF(DFTYPE.EQ.P86)   THEN
         NCORFG    = 9
         EXCHR     = 'HFX'
         CORCHR    = CORNM(9)
         DFTTYP(3) = 1.0D+00
      ELSEIF(DFTYPE.EQ.PZ81)   THEN
         NCORFG    = 10
         EXCHR     = 'HFX'
         CORCHR    = CORNM(10)
         DFTTYP(3) = 1.0D+00
      ELSEIF(DFTYPE.EQ.CTRUHLAR05) THEN
         NCORFG    = 11
         EXCHR     = 'HFX'
         CORCHR    = CORNM(12)
         DFTTYP(3) = 1.0D+00
         IM05      = 1
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.CTRUHLAR052X) THEN
         NCORFG    = 11
         EXCHR     = 'HFX'
         CORCHR    = 'M05-2X'
         DFTTYP(3) = 1.0D+00
         IM05      = 2
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.CVS98) THEN
         NCORFG    = 12
         EXCHR     = 'HFX'
         CORCHR    = CORNM(13)
         DFTTYP(3) = 1.0D+00
         IM06      = 1
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.CTRUHLAR06) THEN
         NCORFG    = 13
         EXCHR     = 'HFX'
         CORCHR    = CORNM(14)
         DFTTYP(3) = 1.0D+00
         IM06      = 3
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.CTRUHLAR06L) THEN
         NCORFG    = 13
         EXCHR     = 'HFX'
         CORCHR    = 'M06-L'
         DFTTYP(3) = 1.0D+00
         IM06      = 1
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.CTRUHLAR062X) THEN
         NCORFG    = 13
         EXCHR     = 'HFX'
         CORCHR    = 'M06-2X'
         DFTTYP(3) = 1.0D+00
         IM06      = 4
         NEEDTAU   = .TRUE.
      ELSEIF(DFTYPE.EQ.CTRUHLAR06HF) THEN
         NCORFG    = 13
         EXCHR     = 'HFX'
         CORCHR    = 'M06-HF'
         DFTTYP(3) = 1.0D+00
         IM06      = 2
         NEEDTAU   = .TRUE.
C
C     ----- EXCHANGE/CORRELATION METHOD -----
C
      ELSEIF(DFTYPE.EQ.PW91) THEN
         NEXFG     = 3
         NCORFG    = 3
         EXCHR     = EXNM(3)
         CORCHR    = CORNM(3)
      ELSEIF(DFTYPE.EQ.PBE) THEN
         NEXFG     = 4
         NCORFG    = 4
         EXCHR     = EXNM(4)
         CORCHR    = CORNM(4)
      ELSEIF(DFTYPE.EQ.TRUHLAR05) THEN
         NEXFG     = 11
         NCORFG    = 11
         EXCHR     = EXNM(8)
         CORCHR    = CORNM(12)
         NEEDTAU   = .TRUE.
         DFTTYP(3) = 0.28D+00
         CM05XF    = ONE-DFTTYP(3)
         IM05      = 1
      ELSEIF(DFTYPE.EQ.TRUHLAR052X) THEN
         NEXFG     = 11
         NCORFG    = 11
         EXCHR     = EXNM(8)
         CORCHR    = CORNM(12)
         NEEDTAU   = .TRUE.
         DFTTYP(3) = 0.56D+00
         CM05XF    = ONE-DFTTYP(3)
         IM05      = 2
      ELSEIF(DFTYPE.EQ.VS98) THEN
         NEXFG     = 12
         NCORFG    = 12
         EXCHR     = EXNM(9)
         CORCHR    = CORNM(13)
         NEEDTAU   = .TRUE.
         IM06      = 1
      ELSEIF(DFTYPE.EQ.TRUHLAR06) THEN
         NEXFG     = 13
         NCORFG    = 13
         EXCHR     = EXNM(10)
         CORCHR    = CORNM(14)
         NEEDTAU   = .TRUE.
         DFTTYP(3) = 0.27D+00
         IM06      = 3
      ELSEIF(DFTYPE.EQ.TRUHLAR06HF) THEN
         NEXFG     = 13
         NCORFG    = 13
         EXCHR     = 'M06-HF'
         CORCHR    = 'M06-HF'
         NEEDTAU   = .TRUE.
         DFTTYP(3) = ONE
         IM06      = 2
      ELSEIF(DFTYPE.EQ.TRUHLAR062X) THEN
         NEXFG     = 13
         NCORFG    = 13
         EXCHR     = 'M06-2X'
         CORCHR    = 'M06-2X'
         NEEDTAU   = .TRUE.
         DFTTYP(3) = 0.54D+00
         IM06      = 4
      ELSEIF(DFTYPE.EQ.TRUHLAR06L) THEN
         NEXFG     = 13
         NCORFG    = 13
         EXCHR     = 'M06-L'
         CORCHR    = 'M06-L'
         NEEDTAU   = .TRUE.
         IM06      = 1
C
C     ----- HYBRID METHOD -----
C
      ELSEIF(DFTYPE.EQ.B3LYP)  THEN
         NEXFG     =-1
         NCORFG    =-1
         EXCHR     = 'B88&HFX'
         CORCHR    = 'LYP88&VWN5'
         CSLT      = 0.08D+00
         CB88      = 0.72D+00
         CLYP      = 0.81D+00
         CVWN      = 0.19D+00
         LCFLAG=.FALSE.
         DFTTYP(3) = 0.20D+00
      ELSEIF(DFTYPE.EQ.B3LYP1) THEN
         NEXFG     =-1
         NCORFG    =-2
         EXCHR     = 'B88&HFX'
         CORCHR    = 'LYP88&VWN1'
         CSLT      = 0.08D+00
         CB88      = 0.72D+00
         CLYP      = 0.81D+00
         CVWN      = 0.19D+00
         DFTTYP(3) = 0.20D+00
       ELSEIF(DFTYPE.EQ.X3LYP) THEN
         NEXFG     =-2
         NCORFG    =-2
         EXCHR     = 'B&PW91&HFX'
         CORCHR    = 'LYP88&VWN1'
         CSLT      = 0.073D+00
         CLYP      = 0.871D+00
         CVWN      = 0.129D+00
         CB88      = 0.542D+00
         CPW91     = 0.167D+00
         DFTTYP(3) = 0.218D+00
C-----------------------------------------------------------------------
C     THERE NEEDS TO BE SOME CLARIFICATION OF THE X3LYP PREFACTORS THAT
C     ARE USED ABOVE.
C     FIRST OFF THE X3LYP PAPER IS:
C     J. CHEM. PHYS. 122 (2005) 014105
C
C     OR IF YOU PREFER PNAS:
C     PNAS 101 (2004) 2673
C
C     IN ANY CASE THE PARAMETERS ARE:
C     PREFACTOR FOR (HF     EXCHANGE) = 0.218 = DFTTYP(3)
C     PREFACTOR FOR (SLATER EXCHANGE) = 0.782 = CSLT
C     PREFACTOR FOR (B88    EXCHANGE) = 0.542 = CB88
C     PREFACTOR FOR (PW91   EXCHANGE) = 0.167 = CPW91
C     PREFACTOR FOR (VWN    CORR.   ) = 0.129 = CVWN
C     PREFACTOR FOR (LYP    CORR.   ) = 0.871 = CLYP
C
C     BUT HERE IS SOMETHING FOR THE NEWBIES.  THE PREFACTOR USED IN THE
C     CODE FOR THE SLATER EXCHANGE IS ACTUALLY:
C
C     CSLT = PREFACTOR FOR (SLATER EXCHANGE) - (CB88 + CPW91)
C
C     JUST TO BEAT A DEAD HORSE SILLY:
C
C     CSLT = 0.782 - (0.542 + 0.129) = 0.073
C-----------------------------------------------------------------------
      ELSEIF(DFTYPE.EQ.HALF) THEN
         NEXFG     = 2
         EXCHR     = 'B88&HFX'
         CORCHR    = 'NONE'
         CB88      = 0.5D+00
         LCFLAG=.FALSE.
         DFTTYP(3) = 0.5D+00
      ELSEIF(DFTYPE.EQ.BHHLYP) THEN
         NEXFG     = 2
         NCORFG    = 2
         EXCHR     = 'B88&HFX'
         CORCHR    = 'LYP88'
         CB88      = 0.5D+00
         LCFLAG=.FALSE.
         DFTTYP(3) = 0.5D+00
      ELSEIF(DFTYPE.EQ.PBE0) THEN
         NEXFG     = 4
         NCORFG    = 4
         EXCHR     = 'PBE&HFX'
         CORCHR    = 'PBE'
         CPBEX     = 0.75D+00
         DFTTYP(3) = 0.25D+00
      ELSEIF(DFTYPE.EQ.PBE0M) THEN
         NEXFG     = 4
         NCORFG    = 3
         EXCHR     = 'PBE&HFX'
         CORCHR    = 'PW91'
         CPBEX     = 0.75D+00
         DFTTYP(3) = 0.25D+00
         CPW91     = ONE
C
C     ----- LC METHOD -----
C
      ELSEIF(DFTYPE.EQ.LCBOP)  THEN
         NEXFG     = 2
         NCORFG    = 7
         EXCHR     = 'B88&LCHFX'
         CORCHR    = CORNM(7)
         QOP       = 2.3670D+00
         LCFLAG=.TRUE.
         DFTTYP(3) = 1.0D+00
      ELSEIF(DFTYPE.EQ.LCBLYP) THEN
         NEXFG     = 2
         NCORFG    = 2
         EXCHR     = 'B88&LCHFX'
         CORCHR    = CORNM(2)
         LCFLAG=.TRUE.
         DFTTYP(3) = 1.0D+00
      ELSE
C
C     ---- OTHER FUNCTIONAL  -----
C     THIS PART GENERATES A ZILLION OR SO COMBINATIONS, MANY OF
C     WHICH ARE SILLY (BUT SOME ARE NOT) BY PUTTING TOGETHER THE
C     SHORT SYMBOLS IN -EX- WITH SHORT SYMBOLS IN -COR-, TO MAKE
C
         CALL NAMEXC(NEX,NCOR,EX,COR,NUMEX,NUMCOR,DFTYPE,
     *               EXNM,CORNM,EXCHR,CORCHR,DQOP,PFTYP)
      END IF
C
      IF(LCFLAG) DFTTYP(3) = 1.0D+00
C
      IF(MASWRK) THEN
         IF(SG1) THEN
            WRITE(IW,9010) DFTYPE,SWOFF,DFTTHR,DFTGTHR
         ELSE
            IF(NLEB(1).EQ.0) THEN
               WRITE (IW,9020) DFTYPE,NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                         SWOFF,SW0,DFTTHR,DFTGTHR
            ELSE
               WRITE (IW,9030) DFTYPE,NRAD,NLEB(1),NRAD0,NLEB0(1),
     *                         SWOFF,SW0,DFTTHR,DFTGTHR
            ENDIF
         ENDIF
         IF(QOP.NE.ZERO) WRITE(IW,9040) QOP,PFTYP
         IF(LCFLAG)      WRITE(IW,9050) LCFLAG
         IF(EXETYP.EQ.CHECK) WRITE(IW,9060)
      ENDIF
      RETURN
C
 9010 FORMAT(/5X,'STANDARD GRID ONE (SG-1) PRUNED GRID OPTIONS'/
     *     5X,43(1H-)/
     *     5X,'DFTTYP=',A8/
     *     5X,'GRID HAS 50 RADIAL POINTS AND AT MOST'/
     *     5X,' 194 ANGULAR POINTS, BEFORE PRUNING'/
     *     5X,'SWOFF =',1P,E12.2,' (PURE SCF -> DFT)'/
     *     5X,'COARSE GRID SWITCHING IS NOT AVAILABLE'/
     *     5X,'THRESH=',1P,E12.2,' GTHRE=',E12.2)
 9020 FORMAT(/5X,'POLAR COORDINATE GRID-BASED DFT OPTIONS'/5X,39(1H-)/
     *     5X,'DFTTYP=',A8/
     *     5X,'NRAD  =',I8,5X,'NTHE  =',I8,5X,'NPHI  =',I8/
     *     5X,'NRAD0 =',I8,5X,'NTHE0 =',I8,5X,'NPHI0 =',I8/
     *     5X,' SWOFF=',1P,E12.2,' (PURE SCF -> DFT)'/
     *     5X,'SWITCH=',1P,E12.2,' (COARSE -> TIGHT GRID)'/
     *     5X,'THRESH=',1P,E12.2,' GTHRE=',E12.2)
 9030 FORMAT(/5X,'LEBEDEV GRID-BASED DFT OPTIONS'/5X,30(1H-)/
     *     5X,'DFTTYP=',A8/
     *     5X,'NRAD  =',I8,5X,'NLEB  =',I8/
     *     5X,'NRAD0 =',I8,5X,'NLEB0 =',I8/
     *     5X,'SWOFF =',1P,E12.2,' (PURE SCF -> DFT)'/
     *     5X,'SWITCH=',1P,E12.2,' (COARSE -> TIGHT GRID)'/
     *     5X,'THRESH=',1P,E12.2,' GTHRE=',E12.2)
 9040 FORMAT(5X,7HQOP   =,F8.3,5X,7HPFTYP =,A8)
 9050 FORMAT(5X,'LONG RANGE CORRECTION=',L2)
 9060 FORMAT(/10X,'**** NOTE - FOR SPEED REASONS,'/
     *       1X,'EXETYP=CHECK RUNS SKIP BECKE ATOMIC PARTITIONING,'/
     *       1X,'LEADING TO RETENTION OF ALL DFT GRID POINTS, AND'/
     *       1X,'THUS AN OVERESTIMATION OF THE ACTUAL MEMORY NEEDS.'//
     *       1X,'MEMORY NEEDS MAY BE ESTIMATED ABOUT 20% TOO HIGH.'//)
C
C     =========================
      ENTRY RNAMEXC(XCHNG,CRRN)
C     =========================
      XCHNG = EXCHR
      CRRN  = CORCHR
      RETURN
      END
C*MODULE DFTXCA  *DECK INPPFREE
      SUBROUTINE INPPFREE(PFTYP)
C**********************************************************************
C     KINFG;
C     1  ...  THOMAS-FERMI-WEIZSACKER (TFW)
C     2  ...  THOMAS-FERMI-LAMDA-WEIZSACKER (TFLW) [LAMBDA = 0.8945]
C     3  ...  LEE-LEE-PARR (LLP)
C     4  ...  THAKKER (THAKKER)
C     5  ...  ??
C     6  ...  PADE (PADE)
C     7  ...  NEW-LLP (NLLP)
C**********************************************************************
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION LLP,NLLP
      SAVE PFKIN
C
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      CHARACTER*10 PFKIN,PFCHR
C
      DATA TFW    /8HTFW     /,TFLW   /8HTFLW    /,LLP    /8HLLP     /
      DATA THAK   /8HTHAK    /,FIT    /8HFIT     /,PADE   /8HPADE    /
      DATA NLLP   /8HNLLP    /
C
      NPFFG=0
C
      IF(PFTYP.EQ.TFW) THEN
         NPFFG  = 1
         PFKIN   = 'TF-W'
         QOP     = 2.3866D+00
      ELSEIF(PFTYP.EQ.TFLW)  THEN
         NPFFG  = 2
         PFKIN   = 'TF-LAMBDAW'
         QOP     = 2.3687D+00
      ELSEIF(PFTYP.EQ.LLP)  THEN
         NPFFG  = 3
         PFKIN   = 'LLP'
         QOP     = 2.3924D+00
      ELSEIF(PFTYP.EQ.THAK)  THEN
         NPFFG  = 4
         PFKIN   = 'THAKKER'
         QOP     = 2.3895D+00
      ELSEIF(PFTYP.EQ.FIT)  THEN
         NPFFG  = 5
         PFKIN   = 'FIT'
         QOP     = 2.3665D+00
      ELSEIF(PFTYP.EQ.PADE)  THEN
         NPFFG  = 6
         PFKIN   = 'PADE'
         QOP     = 2.3833D+00
      ELSEIF(PFTYP.EQ.NLLP)  THEN
         NPFFG  = 7
         PFKIN   = 'NLLP'
         QOP     =2.3653D+00
      ELSE
         WRITE (IW,*) 'PFTYP IS EITHER NOT GIVEN OR INCORRECT'
         CALL ABRT
      ENDIF
      RETURN
C
      ENTRY REFPFREE(PFCHR)
      PFCHR = PFKIN
      RETURN
      END
C*MODULE DFTXCA  *DECK KINETIC
      SUBROUTINE KINETIC(X,X2,TAUS,TAU,DTADX)
C**********************************************************************
C     DETERMINE THE VALUE OF THE TAU FOR EACH KINETIC ENERGY DENSITY
C
C     KINFG;
C     1  ...  THOMAS-FERMI-WEIZSACKER (TFW)
C     2  ...  THOMAS-FERMI-LAMDA-WEIZSACKER (TFLW) [LAMBDA = 0.8945]
C     3  ...  LEE-LEE-PARR (LLP)
C     4  ...  THAKKER (THAKKER)
C     5  ...  ??
C     6  ...  PADE (PADE)
C     7  ...  NEW-LLP (NLLP)
C**********************************************************************
C     THIS ROUTINE MAY BE MISLEADING. NONE OF THESE FUNCTIONALS ARE
C     TRUE META-GGA FUNCTIONALS.
C**********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER ( ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00)
C
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
C     ---- TFW ----
C
      IF (NPFFG.EQ.1) THEN
         TAU  =TAUS+X2/36.0D+00
         DTADX= X/18.0D+00
C
C     ---- TFLW ----
C
      ELSEIF (NPFFG.EQ.2) THEN
         TAU  =TAUS+0.8945D+00*X2/36.0D+00
         DTADX=0.8945D+00*X/18.0D+00
C
C     ----- LLP ----
C
      ELSEIF (NPFFG.EQ.3) THEN
         BB   =0.0044188D+00
         CC   =0.0253D+00
         Q    =SQRT(ONE+X2)
         S    =LOG(X+Q)
         T    =ONE+CC*X*S
         G    =BB*X2/T
         TAU  =TAUS*(ONE+G)
         DTADX=TAUS*(TWO*BB*X+BB*CC*X2*(S-X/Q))/T/T
C
C     ---- THAKKER ----
C
      ELSE IF (NPFFG.EQ.4) THEN
         BB   =0.0055D+00
         CC   =0.0253D+00
         DD   =0.072D+00
         Q    =SQRT(ONE+X2)
         S    =LOG(X+Q)
         T    =ONE+CC*X*S
         G    =BB*X2/T
         H    =ONE+TWO**(5.0D+00/THREE)*X
         TAU  =TAUS*(ONE+G-DD*X/H)
         DTADX=TAUS*(TWO*BB*X+BB*CC*X2*(S-X/Q))/T/T
     >        -TAUS*DD/H/H
C
C     ---- FIT ----
C
      ELSE IF (NPFFG.EQ.4) THEN
         DENOM=ONE+0.0024D+00*X2
         TAU  =TAUS+X2/36.0D+00/DENOM
         DTADX=X/18.0D+00/DENOM-X2*X*0.0048D+00/36.0D+00/DENOM/DENOM
C
C     ---- PADE ----
C
      ELSE IF (NPFFG.EQ.5) THEN
         YYA  = X2/36.0D+00/TAUS
         YYA2 = YYA*YYA
         YYA3 = YYA*YYA2
         YYA4 = YYA2*YYA2
         YYA5 = YYA2*YYA3
         YYA6 = YYA3*YYA3
         A1   = 0.95D+00
         A2   = 14.28111D+00
         A3   =-19.57962D+00
         B1   =-0.05D+00
         B2   = 9.99802D+00
         B3   = 2.96085D+00
         DENOM= B3*YYA3+B2*YYA2+B1*YYA+ONE
         TAU  = TAUS*(9.0D+00*B3*YYA4+A3*YYA3+A2*YYA2+A1*YYA+ONE)/DENOM
         PNOMER= 9.0D+00*B3*B3*YYA6+18.0D+00*B3*B2*YYA5
     >         +(27.0D+00*B3*B1+A3*B2-B3*A2)*YYA4
     >         +(36.0D+00*B3+TWO*(A3*B1-B3*A1))*YYA3
     >         +(THREE*(A3-B3)+A2*B1-B2*A1)*YYA2
     >         +TWO*(A2-B2)*YYA+A1-B1
         DTADX= PNOMER/DENOM/DENOM*X/18.0D+00
      ELSE
         WRITE(IW,*) 'INVALID NUMBER OF NPFFG =',NPFFG
         CALL ABRT
         STOP
      END IF
C
      RETURN
      END
C
C*MODULE DFTXCA  *DECK LENXC
      SUBROUTINE LENXC(CHAR,NUM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*8 CHAR
      CHARACTER*36 ALPHA
      PARAMETER(ALPHA='ABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789')
C
      DO I=1,8
         NUM=I
         IF (INDEX(ALPHA,CHAR(NUM:NUM)).EQ.0) GO TO 15
      ENDDO
C
      NUM = 8
C
 15   CONTINUE
      NUM=NUM-1
      RETURN
      END
C*MODULE DFTXCA  *DECK LSDSLT
      SUBROUTINE LSDSLT(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA,VXCB,XALPHA,
     *                  SIGKA,SIGKB)
C***********************************************************************
C     SLATER(DIRAC) EXCHANGE FUNCTIONAL
C
C     XALF    ....  EXCHANGE ENERGY OF LDA PART
C     B88X    ....  EXCHANGE ENERGY OF GGA PART
C     DUM     ....  EXCHANGE GRADIENT  TO RHO
C     DUMX    ....  EXCHANGE GRADIENT TO X
C     DUMY    ....  EXCHANGE GRADIENT TO Y
C     DUMZ    ....  EXCHANGE GRADIENT TO Z
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00)
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /NLRC  / LCFLAG,LRINT,EMU,EMU2,LRFILE
      LOGICAL LCFLAG,LRINT
C***********************************************************************
      GRDAA =0.0D+00
      GRDBB =0.0D+00
      GRADXA=0.0D+00
      GRADYA=0.0D+00
      GRADZA=0.0D+00
      GRADXB=0.0D+00
      GRADYB=0.0D+00
      GRADZB=0.0D+00
      DKDXA =0.0D+00
      DKDXB =0.0D+00
       IF (RHFGVB) THEN
         RHO13= ROA**(ONE/THREE)
         FACP = -(6.0D+00/PI)**(ONE/THREE)
         FACE = -(THREE/4.0D+00)*(6.0D+00/PI)**(ONE/THREE)
         PDUMA = FACP*RHO13
         VXCB = VXCB
         RHOE = RHO13**4
         XALFA= FTOTWT*RHOE*FACE*TWO
C     ----- FOR OP CORRELATION -----
         SIGKA= -FACE*TWO
C     ------------------------------
         IF (LCFLAG) THEN
         DKDXA=0.0D+00
            CALL LRCFUN(UROHF,ROA,ROB,GRDAA,GRDBB,
     >           GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >           FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >           XALFA,XALFB,B88XA,B88XB,PDUMA,PDUMXA,PDUMYA,
     >           PDUMZA,PDUMB,PDUMXB,PDUMYB,PDUMZB,FACE)
         END IF
         XALPHA = XALPHA + CSLT*XALFA
         VXCA   = VXCA + CSLT*PDUMA
       ELSE IF (UROHF) THEN
         RHO13A=ROA**(ONE/THREE)
         RHO13B=ROB**(ONE/THREE)
         FACP=-(6.0D+00/PI)**(ONE/THREE)
         FACE=-(THREE/4.0D+00)*(6.0D+00/PI)**(ONE/THREE)
         PDUMA = CSLT*FACP*RHO13A
         PDUMB = CSLT*FACP*RHO13B
         RHOEA= RHO13A**4
         RHOEB= RHO13B**4
         XALFA = CSLT*FTOTWT*RHOEA*FACE
         XALFB = CSLT*FTOTWT*RHOEB*FACE
C     ----- FOR OP CORRELATION -----
         SIGKA= -FACE*TWO
         SIGKB= -FACE*TWO
C     ------------------------------
         IF (LCFLAG) THEN
            DKDXA=0.0D+00
            DKDXB=0.0D+00
C     PERHAPS NEED NOT GRDAA,GRDBB=0.0D+00
            CALL LRCFUN(UROHF,ROA,ROB,GRDAA,GRDBB,
     >           GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >           FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >           XALFA,XALFB,B88XA,B88XB,PDUMA,PDUMXA,PDUMYA,
     >           PDUMZA,PDUMB,PDUMXB,PDUMYB,PDUMZB,FACE)
         END IF
         XALPHA = XALPHA + XALFA + XALFB
         VXCA = VXCA+PDUMA
         VXCB = VXCB+PDUMB
       ENDIF
      RETURN
      END
C*MODULE DFTXCA  *DECK LYPCF
      SUBROUTINE LYPCF(RHFGVB,UROHF,NB,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >                 GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                 FTOTWT,ELYP,
     >                 DUMA,DUMXA,DUMYA,DUMZA,DUMB,DUMXB,DUMYB,DUMZB)
C***********************************************************************
C     LYP CORRELATION FUNCTIONAL
C
C     C. LEE, W. YANG AND R. G. PARR
C     PHYSICAL REVIEW B 37, 785-789 (1988)
C
C     -----
C     INPUT
C     -----
C
C     ROA,B   .... SPIN DENSITIES
C     GRAD**  .... DENSITY GRADIENTS
C     FTOTWT  .... VOLUME OF THIS GRID
C
C     ------
C     OUTPUT
C     ------
C
C     ELYP    ....  CORRELATION ENERGY OF LYP
C     DUM***  ....  FIRST DERIVATIVES OF THE EXCHANGE FUNCTINALS
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00)
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
C
       FACE = (TWO**(11.0D+00/THREE))*(THREE/10.0D+00)*
     >        (THREE*PI**2)**(TWO/THREE)
       FACP = ONE*CLYP
       ALYPA=0.04918D+00
       ALYPB=0.13200D+00
       ALYPC=0.25330D+00
       ALYPD=0.34900D+00
C***********************************************************************
       IF (RHFGVB) THEN
C***********************************************************************
            RHOA  =ROA
            RHOB  =ROB
            GAA   =GRDAA
            GBB   =GRDBB
            GAB   =GRDAB
            RHOT  =RHOA+RHOB
            RHO13 =RHOT**(ONE/THREE)
            RHO13A=RHOA**(ONE/THREE)
            RHO13B=RHOB**(ONE/THREE)
            RHO43 =RHO13**4
            RHO83A=RHO13A**8
            RHO83B=RHO13B**8
            XC    =ALYPC/RHO13
            XD    =ALYPD/RHO13
            RX    =XD/(ONE+XD)
            DELTA =XC+RX
            OMEGA =(ALYPA*ALYPB/ALYPD)*RX*EXP(-XC)/RHO43
            DELTAP=(-DELTA+RX*RX)/THREE/RHOT
            OMEGAP=( DELTA-5.0D+00 )/THREE/RHOT
            RA    =RHOA/RHOT
            RB    =RHOB/RHOT
            RAA   =RA*RA
            RBB   =RB*RB
            RAB   =RA*RB
C***********************************************************************
            FAA   =(RAB/9)*(ONE-THREE*DELTA-(DELTA-11)*RA)
            FBB   =(RAB/9)*(ONE-THREE*DELTA-(DELTA-11)*RB)
            FAB   =(RAB/9)*(47-7*DELTA)
            DDAA  =-OMEGA*(FAA-RBB           )
            DDBB  =-OMEGA*(FBB-RAA           )
            DDAB  =-OMEGA*(FAB-    4.0D+00/THREE)
            T1    =(-4.0D+00*ALYPA/ALYPD)*RX*RAB*RHO43
            T2A   =(-FACE*OMEGA*RAB)*RHO83A
            T2B   =(-FACE*OMEGA*RAB)*RHO83B
            T2    =T2A+T2B
            ELYPI =T1+T2+DDAA*GAA+DDBB*GBB+DDAB*GAB
C***********************************************************************
            DFAADA=-(RAB/9)*((THREE+RA)*DELTAP+
     1                             (DELTA-11)*RB/RHOT)
            DFBBDA=-(RAB/9)*((THREE+RB)*DELTAP-
     1                             (DELTA-11)*RB/RHOT)
            DFABDA=-(RAB/9)*7*DELTAP
            D2DAAA= OMEGAP*DDAA-OMEGA*(DFAADA+RBB*(TWO/RHOT))
            D2DABB= OMEGAP*DDBB-OMEGA*(DFBBDA-RAB*(TWO/RHOT))
            D2DAAB= OMEGAP*DDAB-OMEGA* DFABDA
            DT1DA =(T1/RHOT)*(RX/THREE+ONE)
            DT2DA =OMEGAP*(T2A+T2B)+(8.0D+00/THREE)*T2A/RHOA
            DDA   =DT1DA+DT2DA+D2DAAA*GAA
     1                        +D2DABB*GBB
     2                        +D2DAAB*GAB
C***********************************************************************
            ELYP  =ELYP  +FACP*ELYPI*FTOTWT
            U     = DDA
            U     =U*FACP
            DUMA  =DUMA+U
            W     =(DDAA+0.5D+00*DDAB)
            W     =W*FACP
            DUMXA =DUMXA+TWO*W*GRADXA
            DUMYA =DUMYA+TWO*W*GRADYA
            DUMZA =DUMZA+TWO*W*GRADZA
C***********************************************************************
       ELSE IF (UROHF) THEN
C***********************************************************************
            GAA   =GRDAA
            GBB   =GRDBB
            GAB   =GRDAB
            RHOA  =ROA
            RHOB  =ROB
            RHOT  =RHOA+RHOB
            RHO13 =RHOT**(ONE/THREE)
            RHO13A=RHOA**(ONE/THREE)
            RHO13B=RHOB**(ONE/THREE)
            RHO43 =RHO13**4
            RHO83A=RHO13A**8
            RHO83B=RHO13B**8
            XC    =ALYPC/RHO13
            XD    =ALYPD/RHO13
            RX    =XD/(ONE+XD)
            DELTA =XC+RX
            OMEGA =(ALYPA*ALYPB/ALYPD)*RX*EXP(-XC)/RHO43
            DELTAP=(-DELTA+RX*RX)/THREE/RHOT
            OMEGAP=( DELTA-5.0D+00 )/THREE/RHOT
            RA    =RHOA/RHOT
            RB    =RHOB/RHOT
            RAA   =RA*RA
            RBB   =RB*RB
            RAB   =RA*RB
            ZETA  =RA-RB
C***********************************************************************
            FAA   =(RAB/9)*(ONE-THREE*DELTA-(DELTA-11)*RA)
            FBB   =(RAB/9)*(ONE-THREE*DELTA-(DELTA-11)*RB)
            FAB   =(RAB/9)*(47-7*DELTA)
            DDAA  =-OMEGA*(FAA-RBB           )
            DDBB  =-OMEGA*(FBB-RAA           )
            DDAB  =-OMEGA*(FAB-    4.0D+00/THREE)
            T1    =(-4.0D+00*ALYPA/ALYPD)*RX*RAB*RHO43
            T2A   =(-FACE*OMEGA*RAB)*RHO83A
            T2B   =(-FACE*OMEGA*RAB)*RHO83B
            T2    =T2A+T2B
            ELYPI =T1+T2+DDAA*GAA+DDBB*GBB+DDAB*GAB
C***********************************************************************
            DFAADA=-ZETA*FAA/RHOA
     1             -(RAB/9)*((THREE+RA)*DELTAP+
     2                             (DELTA-11)*RB/RHOT)
            DFBBDA=-ZETA*FBB/RHOA
     1             -(RAB/9)*((THREE+RB)*DELTAP-
     2                             (DELTA-11)*RB/RHOT)
            DFABDA=-ZETA*FAB/RHOA
     1             -(RAB/9)*7*DELTAP
            IF(NB.GT.0) THEN
            DFAADB= ZETA*FAA/RHOB
     1             -(RAB/9)*((THREE+RA)*DELTAP-
     2                             (DELTA-11)*RA/RHOT)
            DFBBDB= ZETA*FBB/RHOB
     1             -(RAB/9)*((THREE+RB)*DELTAP+
     2                             (DELTA-11)*RA/RHOT)
            DFABDB= ZETA*FAB/RHOB
     1             -(RAB/9)*7*DELTAP
            ELSE
               DFAADB = ZERO
               DFBBDB = ZERO
               DFABDB = ZERO
            END IF
C***********************************************************************
            D2DAAA= OMEGAP*DDAA-OMEGA*(DFAADA+RBB*(TWO/RHOT))
            D2DBAA= OMEGAP*DDAA-OMEGA*(DFAADB-RAB*(TWO/RHOT))
            D2DABB= OMEGAP*DDBB-OMEGA*(DFBBDA-RAB*(TWO/RHOT))
            D2DBBB= OMEGAP*DDBB-OMEGA*(DFBBDB+RAA*(TWO/RHOT))
            D2DAAB= OMEGAP*DDAB-OMEGA* DFABDA
            D2DBAB= OMEGAP*DDAB-OMEGA* DFABDB
C***********************************************************************
            DT1DA =(T1/RHOT)*(RX/THREE+RB/RA)
            DT2DA =(OMEGAP-ZETA/RHOA)*T2B+
     1             (OMEGAP+(8.0D+00/THREE-ZETA)/RHOA)*T2A
            IF(NB.GT.0) THEN
               DT1DB =(T1/RHOT)*(RX/THREE+RA/RB)
               DT2DB =(OMEGAP+ZETA/RHOB)*T2A+
     1                (OMEGAP+(8.0D+00/THREE+ZETA)/RHOB)*T2B
            ELSE
               DT1DB = ZERO
               DT2DB = ZERO
            END IF
            DDA   =DT1DA+DT2DA+D2DAAA*GAA
     1                        +D2DABB*GBB
     2                        +D2DAAB*GAB
            DDB   =DT1DB+DT2DB+D2DBAA*GAA
     1                        +D2DBBB*GBB
     2                        +D2DBAB*GAB
C***********************************************************************
            ELYP   =ELYP  +FACP*ELYPI*FTOTWT
            UA     = DDA
            UA     =UA*FACP
            DUMA   =DUMA+UA
            UB     = DDB
            UB     =UB*FACP
            DUMB   =DUMB+UB
            WAA    =TWO*DDAA
            WAA    =WAA*FACP
            WBB    =TWO*DDBB
            WBB    =WBB*FACP
            WAB    =    DDAB
            WAB    =WAB*FACP
            DUMXA  =DUMXA+WAA*GRADXA+WAB*GRADXB
            DUMYA  =DUMYA+WAA*GRADYA+WAB*GRADYB
            DUMZA  =DUMZA+WAA*GRADZA+WAB*GRADZB
            DUMXB  =DUMXB+WBB*GRADXB+WAB*GRADXA
            DUMYB  =DUMYB+WBB*GRADYB+WAB*GRADYA
            DUMZB  =DUMZB+WBB*GRADZB+WAB*GRADZA
C***********************************************************************
            ENDIF
        RETURN
      END
C*MODULE DFTXCA  *DECK NAMEXC
      SUBROUTINE NAMEXC(NEX,NCOR,EX,COR,NUMEX,NUMCOR,DFTYPE,
     *                  EXNM,CORNM,EXCHR,CORCHR,DQOP,PFTYP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
C
      CHARACTER*8  EX(NEX),COR(NCOR)
      CHARACTER*10 EXNM(NEX),CORNM(NCOR),EXCHR,CORCHR
      CHARACTER*8  TRYFC,INPXC
C
      DIMENSION NUMEX(NEX),NUMCOR(NCOR),DQOP(NEX)
C
      WRITE(UNIT=INPXC,FMT='(A8)') DFTYPE
C
      DO 10 IEX=1,NEX
         DO 20 ICOR=1,NCOR
            CALL LENXC(EX(IEX),LEX)
            CALL LENXC(COR(ICOR),LCOR)
            TRYFC=EX(IEX)(1:LEX)//COR(ICOR)(1:LCOR)
            IF(TRYFC.EQ.INPXC) GOTO 30
 20      CONTINUE
 10   CONTINUE
C
      WRITE (6,9040) INPXC
 9040 FORMAT(1X,'GRID CODE FINDS UNSUPPORTED DFTTYP=',A8)
      CALL ABRT
      STOP
C
 30   CONTINUE
C          NEXT STATMENT ALLOWS BACKWARDS COMPATIBILITY WITH VWN=VWN5
      IF(ICOR.EQ.6) ICOR=1
      IF(IEX .EQ.7) CSLT = 1.05151D+00
      NEXFG  = NUMEX(IEX)
      NCORFG = NUMCOR(ICOR)
      EXCHR  = EXNM(IEX)
      CORCHR = CORNM(ICOR)
C
      IF (ICOR.EQ.7.AND.QOP.LT.1.0D-15) QOP=DQOP(IEX)
      IF (IEX.EQ.6) CALL INPPFREE(PFTYP)
      RETURN
      END
C*MODULE DFTXCA  *DECK NWPW91CF
      SUBROUTINE NWPW91CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >                    GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                    FTOTWT,ECF,VXCA1,DUMAX,DUMAY,DUMAZ,
     >                    VXCB1,DUMBX,DUMBY,DUMBZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
C
      IF(RHFGVB) THEN
C
         IDERIV   = 1
         NPT      = 1
         RHOA1    = 2.0D+00*ROA
         SIGMAAA1 = 4.0D+00*GRDAA
C
         CALL RKS_PW91C_NWCHEM
     *   (IDERIV,NPT,RHOA1,SIGMAAA1,ZK,VRHOA,VSIGMAAA,V2RHOA2,
     *    V2RHOASIGMAAA,V2SIGMAAA2)
C
         WAA = VSIGMAAA
C
         VXCA1 = VXCA1 + VRHOA
         DUMAX = DUMAX + 1.0D+00*(WAA*GRADXA)
         DUMAY = DUMAY + 1.0D+00*(WAA*GRADYA)
         DUMAZ = DUMAZ + 1.0D+00*(WAA*GRADZA)
C
         ECF = ECF   + ZK * FTOTWT
      END IF
C
      IF(UROHF) THEN
C
         IDERIV   = 1
         NPT      = 1
         RHOA1    = ROA
         SIGMAAA1 = GRDAA
         RHOB1    = ROB
         SIGMABB1 = GRDBB
         SIGMAAB1 = GRDAB
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
              WAA = VSIGMAAA
              WBB = VSIGMABB
              WAB = VSIGMAAB
C
            VXCA1 = VXCA1 + VRHOA
            DUMAX = DUMAX + 2.0D+00*(WAA*GRADXA) + (WAB*GRADXB)*1.0D+00
            DUMAY = DUMAY + 2.0D+00*(WAA*GRADYA) + (WAB*GRADYB)*1.0D+00
            DUMAZ = DUMAZ + 2.0D+00*(WAA*GRADZA) + (WAB*GRADZB)*1.0D+00
C
            VXCB1 = VXCB1 + VRHOB
            DUMBX = DUMBX + 2.0D+00*(WBB*GRADXB) + (WAB*GRADXA)*1.0D+00
            DUMBY = DUMBY + 2.0D+00*(WBB*GRADYB) + (WAB*GRADYA)*1.0D+00
            DUMBZ = DUMBZ + 2.0D+00*(WBB*GRADZB) + (WAB*GRADZA)*1.0D+00
C
              ECF = ECF   + ZK * FTOTWT
C
      END IF
C
      RETURN
      END
C*MODULE DFTXCA  *DECK P86CF
      SUBROUTINE P86CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,VXCB1,DUMBX,DUMBY,DUMBZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
C
C        INTERFACED BY ROBERTO PEVERATI AND KIM BALDRIDGE, U.ZURICH
C
      IF(RHFGVB) THEN
C
         IDERIV   = 1
         NPT      = 1
         RHOA1    = 2.0D+00*ROA
         SIGMAAA1 = 4.0D+00*GRDAA
C
         CALL RKS_C_P86(IDERIV,NPT,RHOA1,SIGMAAA1,ZK,VRHOA,VSIGMAAA,
     &                  V2RHOA2,V2RHOASIGMAAA,V2SIGMAAA2)
C
         WAA = VSIGMAAA
C
         VXCA1 = VXCA1 + VRHOA
         DUMAX = DUMAX + 1.0D+00*(WAA*GRADXA)
         DUMAY = DUMAY + 1.0D+00*(WAA*GRADYA)
         DUMAZ = DUMAZ + 1.0D+00*(WAA*GRADZA)
C
         ECF = ECF   + ZK * FTOTWT
C
      END IF
C
      IF(UROHF) THEN
C
         IDERIV   = 1
         NPT      = 1
         RHOA1    = ROA
         SIGMAAA1 = GRDAA
         RHOB1    = ROB
         SIGMABB1 = GRDBB
         SIGMAAB1 = GRDAB
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
         WAA = VSIGMAAA
         WBB = VSIGMABB
         WAB = VSIGMAAB
C
         VXCA1 = VXCA1 + VRHOA
         DUMAX = DUMAX + 2.0D+00*(WAA*GRADXA) + (WAB*GRADXB)*1.0D+00
         DUMAY = DUMAY + 2.0D+00*(WAA*GRADYA) + (WAB*GRADYB)*1.0D+00
         DUMAZ = DUMAZ + 2.0D+00*(WAA*GRADZA) + (WAB*GRADZB)*1.0D+00
C
         VXCB1 = VXCB1 + VRHOB
         DUMBX = DUMBX + 2.0D+00*(WBB*GRADXB) + (WAB*GRADXA)*1.0D+00
         DUMBY = DUMBY + 2.0D+00*(WBB*GRADYB) + (WAB*GRADYA)*1.0D+00
         DUMBZ = DUMBZ + 2.0D+00*(WBB*GRADZB) + (WAB*GRADZA)*1.0D+00
C
         ECF = ECF   + ZK * FTOTWT
C
      END IF
C
      RETURN
      END
C
C*MODULE DFTXCA  *DECK OPCOR
      SUBROUTINE OPCOR(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                 GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                 FTOTWT,EOP,
     >                 SIGKA,SIGKB,DKDXA,DKDXB,
     >                 DUMA,DUMXA,DUMYA,DUMZA,DUMB,DUMXB,DUMYB,DUMZB)
C***********************************************************************
C     ONE-PARAMETER PROGRESSIVE CORRELATION FUNCTIONAL
C
C     T. TSUNEDA AND K. HIRAO
C     J. CHEM. PHYS. 110, 10664-10678 (1999)
C
C     -----
C     INPUT
C     -----
C
C     ROA,B   .... SPINDENSITIES
C     GRAD**  .... DENSITY GRADIENTS
C     FTOTWT  .... VOLUME OF THIS GRID
C
C     ------
C     OUTPUT
C     ------
C
C     EOP    ....  CORRELATION ENERGY OF OP
C     DUM***  ....  FIRST DERIVATIVES OF THE CORRELATION FUNCTINALS
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
      PARAMETER ( ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00)
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      FACP   = ONE
C***********************************************************************
       IF (RHFGVB) THEN
C***********************************************************************
            RHOI=ROA
            RHOA=ROA
            RHOB=ROA
            RHOP=RHOI**(ONE/THREE)
            RHOE=RHOI*RHOP
            DRHO=SQRT(GRDAA)
            X   =DRHO/RHOE
C***********************************************************************
C      THE EXCHANGE PART DERIVED FROM THE EXCHANGE FUNCTIONAL
C***********************************************************************
            DXDR  =-(4.0D+00/THREE)*X/RHOI
            DKDR  =DKDXA*DXDR
            DKDAA =DKDXA/(TWO*DRHO*RHOE)
C***********************************************************************
C      THE CALCULATION OF THE OP CORRELATION ENERGY
C***********************************************************************
            BETA  =(QOP/TWO)*RHOP*SIGKA
            BETA4 =BETA**4
            BETA3 =BETA**3
            BETA2 =BETA**2
            DENOM =BETA4+(1.1284D+00)*BETA3+(0.3183D+00)*BETA2
            ENEWP =((1.5214D+00)*BETA+(0.5764D+00))/DENOM
            ENEWI =-ENEWP*RHOA*RHOB
            EOP   =EOP  +ENEWI*FTOTWT
C***********************************************************************
C      THE CALCULATION OF THE FIRST DERIVATIVES OF OP
C***********************************************************************
            DEPDB =((-4.5642D+00)*BETA4+(-5.7391D+00)*BETA3
     >             +(-2.4355D+00)*BETA2+(-3.6694D-01)*BETA)/DENOM**2
            DRPKDR=ONE/THREE*RHOI**(-TWO/THREE)*SIGKA+RHOP*DKDR
            DBDR  =QOP/4.0D+00*DRPKDR
            DEDR  =-RHOI*ENEWP-RHOI**2*DBDR*DEPDB
            U     =DEDR
            U     =U*FACP
            DUMA  =DUMA+U
C***********************************************************************
            DBDAA =QOP/4.0D+00*RHOP*DKDAA
            DEDAA =-RHOI**2*DBDAA*DEPDB
            W     =DEDAA
            W     =W*FACP
            DUMXA =DUMXA+TWO*W*GRADXA
            DUMYA =DUMYA+TWO*W*GRADYA
            DUMZA =DUMZA+TWO*W*GRADZA
C***********************************************************************
       ELSE IF (UROHF) THEN
C***********************************************************************
            RHOI  = ROA + ROB
            DRHOA =SQRT(GRDAA)
            DRHOB =SQRT(GRDBB)
            RHOA =ROA
            RHOB =ROB
            IF (RHOA.GT.RHOB*1.0D+04) GOTO 1
            RHO13A=RHOA**(  ONE/THREE)
            RHO13B=RHOB**(  ONE/THREE)
            RHO43A=RHO13A*RHOA
            RHO43B=RHO13B*RHOB
            XA    =DRHOA/RHO43A
            XB    =DRHOB/RHO43B
C***********************************************************************
C      THE EXCHANGE PART DERIVED FROM THE EXCHANGE FUNCTIONAL
C***********************************************************************
            DXADRA  =-(4.0D+00/THREE)*XA/RHOA
            DKADRA  =DKDXA*DXADRA
C
            DXBDRB  =-(4.0D+00/THREE)*XB/RHOB
            DKBDRB  =DKDXB*DXBDRB
C
            DKADGA  =DKDXA/RHO43A
            DKADGAGA=DKADGA/TWO/DRHOA
C
            DKBDGB  =DKDXB/RHO43B
            DKBDGBGB=DKBDGB/TWO/DRHOB
C***********************************************************************
C      THE CALCULATION OF THE OP CORRELATION ENERGY
C***********************************************************************
            BDENOM=RHO13A*SIGKA+RHO13B*SIGKB
            BETA  =QOP*RHO13A*RHO13B*SIGKA*SIGKB/BDENOM
            BETA2 =BETA**2
            BETA3 =BETA**3
            BETA4 =BETA**4
            DENOM =BETA4+(1.1284D+00)*BETA3+(0.3183D+00)*BETA2
            ENEWP =((1.5214D+00)*BETA+0.5764D+00)/DENOM
            ENEWI =-RHOA*RHOB*ENEWP
            EOP   =EOP   +ENEWI *FTOTWT
C***********************************************************************
C      THE CALCULATION OF THE FIRST DERIVATIVES OF OP
C***********************************************************************
            DEPDB  =((-4.5642D+00)*BETA4+(-5.7391D+00)*BETA3
     1              +(-2.4355D+00)*BETA2+(-3.6694D-01)*BETA)
     2              /DENOM**2
            DRPKDRA =ONE/THREE*(RHOA**(-TWO/THREE))*SIGKA
     1               +RHO13A*DKADRA
            DBDRA   =QOP*DRPKDRA*(RHO13B*SIGKB/BDENOM)**2
            DEDRA   =-RHOB*ENEWP-RHOA*RHOB*DBDRA*DEPDB
C***********************************************************************
            DRPKDRB =ONE/THREE*(RHOB**(-TWO/THREE))*SIGKB
     1               +RHO13B*DKBDRB
            DBDRB   =QOP*DRPKDRB*(RHO13A*SIGKA/BDENOM)**2
            DEDRB   =-RHOA*ENEWP-RHOB*RHOA*DBDRB*DEPDB
C***********************************************************************
            DBDGAGA =QOP*RHO13A*(RHO13B**2)*(SIGKB**2)
     1               /(BDENOM**2)*DKADGAGA
            DEDGAGA =-RHOA*RHOB*DBDGAGA*DEPDB
C***********************************************************************
            DBDGBGB =QOP*RHO13B*(RHO13A**2)*(SIGKA**2)
     1           /(BDENOM**2)*DKBDGBGB
            DEDGBGB =-RHOB*RHOA*DBDGBGB*DEPDB
C***********************************************************************
            UA     =DEDRA
            UA     =UA*FACP
            DUMA   =DUMA+UA
            UB     =DEDRB
            UB     =UB*FACP
            DUMB   =DUMB+UB
            WAA    =DEDGAGA
            WAA    =WAA*FACP
            WBB    =DEDGBGB
            WBB    =WBB*FACP
            DUMXA  =DUMXA+TWO*WAA*GRADXA
            DUMYA  =DUMYA+TWO*WAA*GRADYA
            DUMZA  =DUMZA+TWO*WAA*GRADZA
            DUMXB  =DUMXB+TWO*WBB*GRADXB
            DUMYB  =DUMYB+TWO*WBB*GRADYB
            DUMZB  =DUMZB+TWO*WBB*GRADZB
 1          CONTINUE
       ENDIF
C***********************************************************************
      RETURN
      END
C*MODULE DFTXCA  *DECK OPTXXF
      SUBROUTINE OPTXXF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                  GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                  FTOTWT,OPTX,DUMA,DUMXA,DUMYA,DUMZA,
     >                  DUMB,DUMXB,DUMYB,DUMZB)
      IMPLICIT NONE
      COMMON /EXFACT/ COPTX
C***********************************************************************
C     OPTX EXCHANGE FUNCTIONAL OF NICHOLAS C. HANDY AND ARON H. COHEN
C
C     MOL. PHYS. (2001) 99 403-412
C
C     EX[RHO] = -(1.05151 * CX + 1.43169 * U**2) * RHO**(4/3)
C             = -1.05151 * CX * RHO**(4/3) + -1.4369 * U**2 * RHO**(4/3)
C             =  [   LOCAL EXCHANGE      ] + [      GGA EXCHANGE       ]
C             =  [   CACULATED IN LDSLT  ] + [  CALCULATED IN OPTXXF   ]
C
C     CX = DIRAC SLATER COEFFICIENT = (3/4)*(6/PI)**(1/3) IT IS
C          CALCULATED IN THE ROUTINE LSDSLT WHICH WILL EXPLAIN WHY
C          THE VALUE OF -CSLT- IS SET TO 1.05151
C
C     COPTX = CONTAINS THE FACTOR 1.43169
C***********************************************************************
      DOUBLE PRECISION ROA, ROB, GRDAA, GRDBB
      DOUBLE PRECISION GRADXA, GRADYA, GRADZA, GRADXB, GRADYB, GRADZB
      DOUBLE PRECISION DUMA, DUMXA, DUMYA, DUMZA
      DOUBLE PRECISION DUMB, DUMXB, DUMYB, DUMZB
      DOUBLE PRECISION FTOTWT, OPTX
      DOUBLE PRECISION COPTX
C
      LOGICAL RHFGVB, UROHF
C
      DOUBLE PRECISION ONE, TWO, THREE, FOUR, GAMMA
      PARAMETER (ONE=1.0D+00,TWO=2.0D+00,THREE=3.0D+00,FOUR=4.0D+00)
      PARAMETER (GAMMA=0.006D+00)
C
      DOUBLE PRECISION RHO, GRD, RHO13, RHO43, X
      DOUBLE PRECISION U, DUDX, E, DEDRHO, DEDGRD
C
      IF (RHFGVB) THEN
        RHO = ROA
        GRD = GRDAA
        RHO13 = RHO**(ONE/THREE)
        RHO43 = RHO**(FOUR/THREE)
        X = SQRT(GRD)/RHO43
        U = (GAMMA*X*X)/(ONE+(GAMMA*X*X))
        DUDX = (TWO*GAMMA*X)/(ONE+(GAMMA*X*X))**TWO
C
        E = -U*U*RHO43
C       THE EXCHANGE ENERGY IS DEFINED ABOVE IN -E-
C
        OPTX = OPTX+E*TWO*FTOTWT*COPTX
C       NOTICED THAT YOU MUST MULTIPLE THE EXCHANGE ENERGY BY SEVERAL
C       FACTORS:
C          TWO : TAKES INTO ACCOUNT THAT IT IS A CLOSED SHELL SYSTEM
C       FTOTWT : TAKES INTO ACCOUNT THE CURRENT GRID WEIGHT
C        COPTX : TAKES INTO ACCOUNT THE FACTOR 1.43169 IN THE DEFINITION
C
        DEDRHO = (FOUR/THREE)*RHO13*((E/RHO43)-(X*(-TWO*U*DUDX)))
C       THE DERIVATIVE OF THE EXCHANGE ENERGY WITH RESPECT TO THE
C       DENSITY IS DEFINED ABOVE IN -DEDRHO-
C
        DEDGRD = (ONE/TWO)*(-TWO*U*DUDX)/SQRT(GRD)
C       THE DERIVATIVE OF THE EXCHANGE ENERGY WITH RESPECT TO THE
C       DENSITY GRADIENT IS DEFINED ABOVE IN -DEDGRD-
C
        DUMA = DUMA+DEDRHO*COPTX
C       NOTICE THAT THERE IS NO FACTOR OF TWO NOR IS THERE AN FTOTWT
C
        DUMXA = DUMXA+(DEDGRD*GRADXA*COPTX*TWO)
        DUMYA = DUMYA+(DEDGRD*GRADYA*COPTX*TWO)
        DUMZA = DUMZA+(DEDGRD*GRADZA*COPTX*TWO)
C       NOTICE THAT THERE IS NO FTOTWT
C
      ELSE IF (UROHF) THEN
C       THE SPIN UNRESTRICTED CODE JUST CONTAINS TWO CONSECTUTIVE
C       EXECUTION OF THE SPIN RESTRICTED CODE
C
C       ALPHA SPIN
C
        RHO = ROA
        GRD = GRDAA
        RHO13 = RHO**(ONE/THREE)
        RHO43 = RHO**(FOUR/THREE)
        X = SQRT(GRD)/RHO43
        U = (GAMMA*X*X)/(ONE+(GAMMA*X*X))
        DUDX = (TWO*GAMMA*X)/(ONE+(GAMMA*X*X))**TWO
        E = -U*U*RHO43
C
        OPTX = OPTX+E*FTOTWT*COPTX
C       NOTICE THAT THE FACTOR OF TWO DOES NOT APPEAR IN THE SPIN
C       UNRESTRICTED CASE
C
        DEDRHO = (FOUR/THREE)*RHO13*((E/RHO43)-(X*(-TWO*U*DUDX)))
        DEDGRD = (ONE/TWO)*(-TWO*U*DUDX)/SQRT(GRD)
        DUMA = DUMA+DEDRHO*COPTX
        DUMXA = DUMXA+(DEDGRD*GRADXA*COPTX*TWO)
        DUMYA = DUMYA+(DEDGRD*GRADYA*COPTX*TWO)
        DUMZA = DUMZA+(DEDGRD*GRADZA*COPTX*TWO)
C
C       BETA SPIN
C
        RHO = ROB
        GRD = GRDBB
        RHO13 = RHO**(ONE/THREE)
        RHO43 = RHO**(FOUR/THREE)
        X = SQRT(GRD)/RHO43
        U = (GAMMA*X*X)/(ONE+(GAMMA*X*X))
        DUDX = (TWO*GAMMA*X)/(ONE+(GAMMA*X*X))**TWO
        E = -U*U*RHO43
C
        OPTX = OPTX+E*FTOTWT*COPTX
C       NOTICE THAT THE FACTOR OF TWO DOES NOT APPEAR IN THE SPIN
C       UNRESTRICTED CASE
C
        DEDRHO = (FOUR/THREE)*RHO13*((E/RHO43)-(X*(-TWO*U*DUDX)))
        DEDGRD = (ONE/TWO)*(-TWO*U*DUDX)/SQRT(GRD)
        DUMB = DUMB+DEDRHO*COPTX
        DUMXB = DUMXB+(DEDGRD*GRADXB*COPTX*TWO)
        DUMYB = DUMYB+(DEDGRD*GRADYB*COPTX*TWO)
        DUMZB = DUMZB+(DEDGRD*GRADZB*COPTX*TWO)
      END IF
C
      END
C*MODULE DFTXCA  *DECK PBECF
      SUBROUTINE PBECF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,VXCB1,DUMBX,DUMBY,DUMBZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
C
C     THE FOLLOWING IS THE BASIC SETUP FOR INTEGRATING A GGA TYPE
C     CORRELATION FUNCTIONAL FROM THE DFT REPOSITORY.
C
      IF(RHFGVB) THEN
C
         IDERIV   = 1
         NPT      = 1
         RHOA1    = 2.0D+00*ROA
         SIGMAAA1 = 4.0D+00*GRDAA
C
         CALL RKS_C_PBE(IDERIV,NPT,RHOA1,SIGMAAA1,ZK,VRHOA,
     *                  VSIGMAAA,V2RHOA2,V2RHOASIGMAAA,V2SIGMAAA2)
C
         WAA = VSIGMAAA
C
         VXCA1 = VXCA1 + VRHOA
         DUMAX = DUMAX + 1.0D+00*(WAA*GRADXA)
         DUMAY = DUMAY + 1.0D+00*(WAA*GRADYA)
         DUMAZ = DUMAZ + 1.0D+00*(WAA*GRADZA)
C
         ECF = ECF   + ZK * FTOTWT
      END IF
C
      IF(UROHF) THEN
C
         IDERIV   = 1
         NPT      = 1
         RHOA1    = ROA
         SIGMAAA1 = GRDAA
         RHOB1    = ROB
         SIGMABB1 = GRDBB
         SIGMAAB1 = GRDAB
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
              WAA = VSIGMAAA
              WBB = VSIGMABB
              WAB = VSIGMAAB
C
            VXCA1 = VXCA1 + VRHOA
            DUMAX = DUMAX + 2.0D+00*(WAA*GRADXA) + (WAB*GRADXB)*1.0D+00
            DUMAY = DUMAY + 2.0D+00*(WAA*GRADYA) + (WAB*GRADYB)*1.0D+00
            DUMAZ = DUMAZ + 2.0D+00*(WAA*GRADZA) + (WAB*GRADZB)*1.0D+00
C
            VXCB1 = VXCB1 + VRHOB
            DUMBX = DUMBX + 2.0D+00*(WBB*GRADXB) + (WAB*GRADXA)*1.0D+00
            DUMBY = DUMBY + 2.0D+00*(WBB*GRADYB) + (WAB*GRADYA)*1.0D+00
            DUMBZ = DUMBZ + 2.0D+00*(WBB*GRADZB) + (WAB*GRADZA)*1.0D+00
C
              ECF = ECF   + ZK * FTOTWT
C
      END IF
C
      RETURN
      END
C
C    DEPOSITORY'S CODE NOT USED, IN FAVOR OF THE FIXED ORIGINAL ROUTINE
C
C---C*MODULE DFTXCA  *DECK PBEXF
C---      SUBROUTINE PBEXF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
C---     >                 GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
C---     >                 FTOTWT,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
C---     >                 VXCB1,DUMBX,DUMBY,DUMBZ)
C---      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C---      LOGICAL UROHF, RHFGVB
C---      COMMON /CPBE0 / CPBEX
C---C
C---C     THE FOLLOWING IS THE BASIC SETUP FOR INTEGRATING A GGA TYPE
C---C     EXCHANGE FUNCTIONAL FROM THE DFT REPOSITORY.
C---C
C---      IF(RHFGVB) THEN
C---C
C---         IDERIV   = 1
C---         NPT      = 1
C---         RHOA1    = 2.0D+00*ROA
C---         SIGMAAA1 = 4.0D+00*GRDAA
C---C
C---        CALL RKS_X_PBE
C---     & (IDERIV,NPT,RHOA1,SIGMAAA1,
C---     &  ZK,VRHOA,VSIGMAAA,
C---     &  V2RHOA2,V2RHOASIGMAAA,V2SIGMAAA2)
C---C
C---              WAA = VSIGMAAA
C---C
C---            VXCA1 = VXCA1 + VRHOA * CPBEX
C---            DUMAX = DUMAX + 1.0D+00*(WAA*GRADXA) * CPBEX
C---            DUMAY = DUMAY + 1.0D+00*(WAA*GRADYA) * CPBEX
C---            DUMAZ = DUMAZ + 1.0D+00*(WAA*GRADZA) * CPBEX
C---C
C---            XGRD = XGRD + ZK * FTOTWT * CPBEX
C---C
C---      END IF
C---C
C---      IF(UROHF) THEN
C---C
C---         IDERIV   = 1
C---         NPT      = 1
C---         RHOA1    = ROA
C---         SIGMAAA1 = GRDAA
C---         RHOB1    = ROB
C---         SIGMABB1 = GRDBB
C---C
C---        CALL UKS_X_PBE
C---     & (IDERIV,NPT,RHOA1,RHOB1,SIGMAAA1,SIGMABB1,
C---     &  ZK,VRHOA,VRHOB,VSIGMAAA,VSIGMABB,VSIGMAAB,
C---     &  V2RHOA2,V2RHOB2,V2RHOAB,
C---     &  V2RHOASIGMAAA,V2RHOASIGMAAB,V2RHOASIGMABB,
C---     &  V2RHOBSIGMABB,V2RHOBSIGMAAB,V2RHOBSIGMAAA,
C---     &  V2SIGMAAA2,V2SIGMAAAAB,V2SIGMAAABB,
C---     &  V2SIGMAAB2,V2SIGMAABBB,V2SIGMABB2)
C---C
C---              WAA = VSIGMAAA
C---              WBB = VSIGMABB
C---C
C---            VXCA1 = VXCA1 + VRHOA * CPBEX
C---            DUMAX = DUMAX + 2.0D+00*(WAA*GRADXA) * CPBEX
C---            DUMAY = DUMAY + 2.0D+00*(WAA*GRADYA) * CPBEX
C---            DUMAZ = DUMAZ + 2.0D+00*(WAA*GRADZA) * CPBEX
C---C
C---            VXCB1 = VXCB1 + VRHOB * CPBEX
C---            DUMBX = DUMBX + 2.0D+00*(WBB*GRADXB) * CPBEX
C---            DUMBY = DUMBY + 2.0D+00*(WBB*GRADYB) * CPBEX
C---            DUMBZ = DUMBZ + 2.0D+00*(WBB*GRADZB) * CPBEX
C---C
C---            XGRD = XGRD + ZK * FTOTWT * CPBEX
C---C
C---      END IF
C---C
C---      END
C
C   THIS ROUTINE IS KEPT BECAUSE IT HAS PARTS RELATING TO THE OP
C   CORRELATION FUNCTIONAL WHICH ARE LOST IF USING THE REPOSITORY CODE.
C   NOTE THAT THIS HAS HAD A CONSTANT CHANGED (CMU) TO FIX SOME
C   SMALL NUMERICAL DISCREPANCIES FOUND IN THE ORIGINAL CODE.
C
C*MODULE DFTXCA  *DECK PBEXFOP
      SUBROUTINE PBEXFOP(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                   GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                   FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >                   XALF,B88X,DUMA,DUMXA,DUMYA,DUMZA,
     >                   DUMB,DUMXB,DUMYB,DUMZB)
C******************************************************************
C
C     PERDEW-BURKE-ERNZERHOF (PBE) EXCHANGE FUNCTIONAL
C
C     J. P. PERDEW , K. BURKE AND M. ERNZERHOF
C     PHYS. REV. LET. 77,3865 (1996)
C
C     -----
C     INPUT
C     -----
C
C     ROA,B   .... SPINDENSITIES
C     GRAD**  .... DENSITY GRADIENTS
C     FTOTWT  .... VOLUME OF THIS GRID
C
C     ------
C     OUTPUT
C     ------
C
C     XALF    ....  EXCHANGE ENERGY OF LDA PART
C     B88X    ....  EXCHANGE ENERGY OF GGA PART
C     DUM***  ....  FIRST DERIVATIVES OF THE EXCHANGE FUNCTINALS
C     SIGKA,B ....  EXC=-1/2*INT(SIGK*RHO^(4/3))
C     DKDXA,B ....  DSIGK/DX
C
C******************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
      COMMON /CPBE0 / CPBEX
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
C******************************************************************
       X1 = (3.0D+00/PI)**(ONE/3.0D+00)
       X2 = (3.0D+00*PI*PI)**(ONE/3.0D+00)
       TWO13 = TWO**(ONE/3.0D+00)
       FACE  =-0.75D+00*X1*TWO13
       FACP  = 4.0D+00/3.0D+00
       CMU   = 0.2195149727645171D+00
       IF (RHFGVB) THEN
         RHO=ROA
         GRD=GRDAA
         RHO13=RHO**(ONE/3.0D+00)
         RHOE =RHO13**4
         DRHO =SQRT(GRD)
         X    =DRHO/RHOE
         S    =X/TWO/TWO13/X2
         Q    =ONE+CMU*S*S/0.804D+00
         G    =FACE*0.804D+00*(ONE-ONE/Q)
         H    =FACE*TWO*CMU*S/Q**2
         R    =FACP*RHO13
         E    =FACE
         XCL  =E*RHOE
         XCNL =G*RHOE
         XALF =XALF +CB88*XCL *TWO*FTOTWT*CPBEX
         B88X =B88X +CB88*XCNL*TWO*FTOTWT*CPBEX
         U    =CB88*R*(E+G-S*H)
         DUMA =DUMA+U*CPBEX
         W    =CB88*0.5D+00*H/TWO/TWO13/X2/SQRT(GRD)
         DUMXA=DUMXA+TWO*W*GRADXA*CPBEX
         DUMYA=DUMYA+TWO*W*GRADYA*CPBEX
         DUMZA=DUMZA+TWO*W*GRADZA*CPBEX
C     ----- FOR OP CORRELATION -----
         SIGKA= -(FACE+G)*TWO
         DKDXA= -H/TWO13/X2
C     -----------------------------
C******************************************************************
C      OPEN SHELL CALCULATION
C******************************************************************
       ELSE IF (UROHF) THEN
         RHO13A=ROA**(ONE/3.0D+00)
         RHOEA =RHO13A**4
         SA    =SQRT(GRDAA)/RHOEA/TWO/TWO13/X2
         QA    =ONE+CMU*SA*SA/0.804D+00
         GA    =FACE*0.804D+00*(ONE-ONE/QA)
         HA    =FACE*TWO*CMU*SA/QA**2
         RA    =FACP*RHO13A
         EA    =FACE
         UA    =CSLT*RA*(EA+GA-SA*HA)
         DUMA  =DUMA+UA*CPBEX
         WA    =CB88*0.5D+00*HA/TWO/TWO13/X2/SQRT(GRDAA)
         DUMXA =DUMXA+TWO*WA*GRADXA*CPBEX
         DUMYA =DUMYA+TWO*WA*GRADYA*CPBEX
         DUMZA =DUMZA+TWO*WA*GRADZA*CPBEX
         RHO13B=ROB**(ONE/3.0D+00)
         RHOEB =RHO13B**4
         EB    =FACE
         IF(ROB.LT.1.0D-15) THEN
           GB  = ZERO
           GOTO 9999
         ENDIF
         SB    =SQRT(GRDBB)/RHOEB/TWO/TWO13/X2
         QB    =ONE+CMU*SB*SB/0.804D+00
         GB    =FACE*0.804D+00*(ONE-ONE/QB)
         HB    =FACE*TWO*CMU*SB/QB**2
         RB    =FACP*RHO13B
         UB    =CSLT*RB*(EB+GB-SB*HB)
         DUMB  =DUMB+UB*CPBEX
         WB    =CB88*0.5D+00*HB/TWO/TWO13/X2/SQRT(GRDBB)
         DUMXB =DUMXB+TWO*WB*GRADXB*CPBEX
         DUMYB =DUMYB+TWO*WB*GRADYB*CPBEX
         DUMZB =DUMZB+TWO*WB*GRADZB*CPBEX
 9999    CONTINUE
         XCL   =EA*RHOEA+EB*RHOEB
         XCNL  =GA*RHOEA+GB*RHOEB
         XALF  =XALF  +CSLT*XCL *FTOTWT*CPBEX
         B88X  =B88X  +CB88*XCNL*FTOTWT*CPBEX
C      ----- FOR OP CORRELATION -----
         SIGKA= -(FACE+GA)*TWO
         SIGKB= -(FACE+GB)*TWO
         DKDXA= -HA/TWO13/X2
         DKDXB= -HB/TWO13/X2
C     -------------------------------
       ENDIF
C*******************************************************************
      RETURN
      END
C
C*MODULE DFTXCA  *DECK PFREE
      SUBROUTINE PFREE(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >                 GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                 FTOTWT,SIGKA,SIGKB,DKDXA,DKDXB,
     >                 XALF,PFX,DUMA,DUMXA,DUMYA,DUMZA,
     >                 DUMB,DUMXB,DUMYB,DUMZB)
C
C***********************************************************************
C
C     PARAMETER-FREE EXCHANGE FUNCTIONAL
C
C     T. TSUNEDA AND K. HIRAO
C     PHYS. REV. B IN PRESS (2000)
C
C     -----
C     INPUT
C     -----
C
C     ROA,B   .... SPIN DENSITIES
C     GRAD**  .... DENSITY GRADIENTS
C     TOTWT   .... VOLUME OF THIS GRID
C
C     ------
C     OUTPUT
C     ------
C
C     PFX    ....  EXCHANGE ENERGY OF GGA PART
C     DUM***  ....  FIRST DERIVATIVES OF THE EXCHANGE FUNCTINALS
C     SIGKA,B ....  EXC=-1/2*INT(SIGK*RHO^(4/3))
C     DKDXA,B ....  DSIGK/DX
C
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (ZERO=0.0D+00,HALF=0.5D+00,
     >           ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00)
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      LOGICAL UROHF, RHFGVB
C
C     ---- CALCULATE PHYSICAL CONSTANTS ----
C
      FACP  = 4.0D+00/3.0D+00
      X1    = (3.0D+00/PI)**(ONE/3.0D+00)
      FACE  = 1.5D+00*X1*(ONE/TWO )**(TWO/THREE)
      TAUS  =(THREE/5.0D+00)*(6.0D+00*PI*PI)**(TWO/THREE)
      COEF1 = 27.0D+00*PI/5.0D+00
      COEF2 = 7.0D+00/108.0D+00
C
C***********************************************************************
C      CLOSED SHELL CALCULATION
C***********************************************************************
      IF (RHFGVB) THEN
         RHO = ROA
         GRD = GRDAA
         RHO13 = RHO**(ONE/THREE)
         RHOE = RHO*RHO13
         X = SQRT(GRD)/RHOE
         X2 = X*X
C
C     ---- CALCULATE KINETIC DENSITY -----
C
         CALL KINETIC(X,X2,TAUS,TAU,DTADX)
C
C     ---- CALCULATE K-SIGMA AND IT'S DERIVATIVE -----
C
         SIGKA =-HALF*COEF1/TAU*(ONE+COEF2*X2/TAU)
         DKDXA = HALF*COEF1/TAU*(DTADX+TWO*COEF2*X*(X*DTADX/TAU-ONE))
     >        /TAU
C
         R    =FACP*RHO13
         XCNL =SIGKA*RHOE
         PFX  =PFX +CB88*TWO*XCNL*FTOTWT
         U    =CB88*R*(SIGKA-X*DKDXA)
         DUMA =DUMA+U
         W    =CB88*HALF*DKDXA/SQRT(GRD)
         DUMXA=DUMXA+TWO*W*GRADXA
         DUMYA=DUMYA+TWO*W*GRADYA
         DUMZA=DUMZA+TWO*W*GRADZA
         XALF =ZERO
C     ---- FOR OP CORRELATION ----
         IF (-SIGKA.LE.FACE) THEN
            SIGKA = -FACE
            DKDXA = ZERO
         ENDIF
         SIGKA=-SIGKA-SIGKA
         DKDXA=-DKDXA-DKDXA
C
C***********************************************************************
C     OPEN SHELL CALCULATION
C***********************************************************************
      ELSE IF (UROHF) THEN
         RHO13A = ROA**(ONE/THREE)
         RHOEA  = ROA*RHO13A
         XA     = SQRT(GRDAA)/RHOEA
         X2A    = XA*XA
C
C     ---- CALCULATE KINETIC DENSITY -----
C
         CALL KINETIC(XA,X2A,TAUS,TAUA,DTADXA)
C
C     ---- CALCULATE K-SIGMA AND IT'S DERIVATIVE -----
C
         SIGKA =-HALF*COEF1/TAUA*(ONE+COEF2*X2A/TAUA)
         DKDXA = HALF*COEF1/TAUA*(DTADXA+TWO*COEF2*XA
     >        *(XA*DTADXA/TAUA-ONE))/TAUA
C
         RA    = FACP*RHO13A
         UA    = CB88*RA*(SIGKA-XA*DKDXA)
         DUMA  = DUMA+UA
         WA    = CB88*HALF*DKDXA/SQRT(GRDAA)
         DUMXA = DUMXA+TWO*WA*GRADXA
         DUMYA = DUMYA+TWO*WA*GRADYA
         DUMZA = DUMZA+TWO*WA*GRADZA
C
C     ---- BETA CONTRIBUTION
C
         IF(ROB.LT.1.0D-15) THEN
            SIGKB=ZERO
           GOTO 9999
         ENDIF
C
         RHO13B = ROB**(ONE/THREE)
         RHOEB  = ROB*RHO13B
         XB     = SQRT(GRDBB)/RHOEB
         X2B    = XB*XB
C
C     ---- CALCULATE KINETIC DENSITY -----
C
         CALL KINETIC(XB,X2B,TAUS,TAUB,DTADXB)
C
C     ---- CALCULATE K-SIGMA AND IT'S DERIVATIVE -----
C
         SIGKB =-HALF*COEF1/TAUB*(ONE+COEF2*X2B/TAUB)
         DKDXB = HALF*COEF1/TAUB*(DTADXB+TWO*COEF2*XB
     >        *(XB*DTADXB/TAUB-ONE))/TAUB
C
         RB    = FACP*RHO13B
         UB    = CB88*RB*(SIGKB-XB*DKDXB)
         DUMB  = DUMB+UB
         WB    = CB88*HALF*DKDXB/SQRT(GRDBB)
         DUMXB = DUMXB+TWO*WB*GRADXB
         DUMYB = DUMYB+TWO*WB*GRADYB
         DUMZB = DUMZB+TWO*WB*GRADZB
 9999    CONTINUE
         XCNL  = SIGKA*RHOEA+SIGKB*RHOEB
         PFX   = PFX +CB88*XCNL*FTOTWT
         XALF =ZERO
C     ---- FOR OP CORRELATION ----
         IF (-SIGKA.LE.FACE) THEN
            SIGKA = -FACE
            DKDXA = ZERO
         ENDIF
         IF (-SIGKB.LE.FACE) THEN
            SIGKB = -FACE
            DKDXB = ZERO
         ENDIF
         SIGKA=-SIGKA-SIGKA
         DKDXA=-DKDXA-DKDXA
         SIGKB=-SIGKB-SIGKB
         DKDXB=-DKDXB-DKDXB
C
      END IF
      RETURN
      END
C
C*MODULE DFTXCA  *DECK PW91CF
      SUBROUTINE PW91CF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,GRDAB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,ECF,
     >        VXCA1,DUMAX,DUMAY,DUMAZ,VXCB1,DUMBX,DUMBY,DUMBZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
C
C     THE PW91 CORRELATION FUNCTIONAL.
C
C     PLEASE NOTE, THE DFT REPOSITORY USED THE CONSTANTS
C       EPS-C(R,0)=0.031091     EPS-C(R,1)=0.015545     -AC=0.016887
C     WHEREAS OTHER PROGRAMS, NOTABLY NWCHEM, USE
C                  0.030090700             0.015545350      0.016886900
C     FOR THE REPRESENTATION OF THE ELECTRON-GAS CORRELATION
C     FUNCTIONAL INCLUDED IN THE TOTAL PW91 GGA.  THE CONSTANTS
C     USED HERE ARE ALL OF THE DIGITS APPEARING IN TABLE 1 OF
C     THE RELEVANT DEFINITION:
C         J.P.PERDEW, Y.WANG, PHYS.REV.B 45, 13244-13249(1992)
C
C     OF COURSE, THE SMALL DIFFERENCES IN THE CONSTANTS MAKE A SMALL
C     IMPACT ON THE ENERGIES OF MOLECULES.  FOR EXAMPLE, IN ETHANOL,
C     WHICH HAS THREE NON-HYDROGEN ATOMS, THERE IS A DIFFERENCE
C     BETWEEN THE TWO PROGRAMS OF ABOUT 5D-6 HARTREE, ONCE VERY
C     FINE GRIDS ARE USED, ALONG WITH GREAT CARE TO REMOVE ALL OTHER
C     DIFFERENCES (GEOMETRY, BASIS, ET CETERA).
C
C     IMPLEMENTING THE PUBLISHED FUNCTIONAL SHOULD USE ANY CONSTANTS,
C     EXACTLY AS PUBLISHED.
C
C     THE FOLLOWING IS THE BASIC SETUP FOR INTEGRATING A GGA TYPE
C     CORRELATION FUNCTIONAL FROM THE DFT REPOSITORY.
C
      IF(RHFGVB) THEN
C
         IDERIV   = 1
         NPT      = 1
         RHOA1    = 2.0D+00*ROA
         SIGMAAA1 = 4.0D+00*GRDAA
C
         CALL RKS_C_PW91(IDERIV,NPT,RHOA1,SIGMAAA1,ZK,VRHOA,
     *                   VSIGMAAA,V2RHOA2,V2RHOASIGMAAA,V2SIGMAAA2)
C
         WAA = VSIGMAAA
C
         VXCA1 = VXCA1 + VRHOA
         DUMAX = DUMAX + 1.0D+00*(WAA*GRADXA)
         DUMAY = DUMAY + 1.0D+00*(WAA*GRADYA)
         DUMAZ = DUMAZ + 1.0D+00*(WAA*GRADZA)
C
         ECF = ECF   + ZK * FTOTWT
      END IF
C
      IF(UROHF) THEN
C
         IDERIV   = 1
         NPT      = 1
         RHOA1    = ROA
         SIGMAAA1 = GRDAA
         RHOB1    = ROB
         SIGMABB1 = GRDBB
         SIGMAAB1 = GRDAB
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
              WAA = VSIGMAAA
              WBB = VSIGMABB
              WAB = VSIGMAAB
C
            VXCA1 = VXCA1 + VRHOA
            DUMAX = DUMAX + 2.0D+00*(WAA*GRADXA) + (WAB*GRADXB)*1.0D+00
            DUMAY = DUMAY + 2.0D+00*(WAA*GRADYA) + (WAB*GRADYB)*1.0D+00
            DUMAZ = DUMAZ + 2.0D+00*(WAA*GRADZA) + (WAB*GRADZB)*1.0D+00
C
            VXCB1 = VXCB1 + VRHOB
            DUMBX = DUMBX + 2.0D+00*(WBB*GRADXB) + (WAB*GRADXA)*1.0D+00
            DUMBY = DUMBY + 2.0D+00*(WBB*GRADYB) + (WAB*GRADYA)*1.0D+00
            DUMBZ = DUMBZ + 2.0D+00*(WBB*GRADZB) + (WAB*GRADZA)*1.0D+00
C
              ECF = ECF   + ZK * FTOTWT
C
      END IF
C
      RETURN
      END
C*MODULE DFTXCA  *DECK PW91LDACF
      SUBROUTINE PW91LDACF(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA1,VXCB1,
     +                     EPW91LDA)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
      COMMON /CPW91L/ CPW91LDA
C
C     THE FOLLOWING IS THE BASIC SETUP FOR INTEGRATING AN LDA TYPE
C     CORRELATION FUNCTIONAL FROM THE DFT REPOSITORY.
C
      IDERIV= 1
      NPT   = 1
C
      IF(RHFGVB)THEN
C
         RHOA1 = ROA + ROB
C
         CALL RKS_C_PW92(IDERIV,NPT,RHOA1,ZK,VRHOA,V2RHOA2)
C
         VXCA1     = VXCA1     +     VRHOA*CPW91LDA
         EPW91LDA  = EPW91LDA  + ZK*FTOTWT*CPW91LDA
C
      END IF
C
      IF(UROHF) THEN
C
         RHOA1 = ROA
         RHOB1 = ROB
C
         CALL UKS_C_PW92(IDERIV,NPT,RHOA1,RHOB1,ZK,
     &                   VRHOA,VRHOB,V2RHOA2,V2RHOB2,V2RHOAB)
C
         VXCA1     = VXCA1     +     VRHOA*CPW91LDA
         VXCB1     = VXCB1     +     VRHOB*CPW91LDA
         EPW91LDA  = EPW91LDA  + ZK*FTOTWT*CPW91LDA
C
      END IF
C
      RETURN
      END
C*MODULE DFTXCA  *DECK PW91XF
      SUBROUTINE PW91XF(RHFGVB,UROHF,ROA,ROB,GRDAA,GRDBB,
     >        GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >        FTOTWT,XPW91,VXCA1,DUMAX,DUMAY,DUMAZ,
     >        VXCB1,DUMBX,DUMBY,DUMBZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
      COMMON /CX3LYP/ CPW91
C
C     THE FOLLOWING IS THE BASIC SETUP FOR INTEGRATING A GGA TYPE
C     EXCHANGE FUNCTIONAL FROM THE DFT REPOSITORY.
C
      IF(RHFGVB) THEN
C
         IDERIV   = 1
         NPT      = 1
         RHOA1    = 2.0D+00*ROA
         SIGMAAA1 = 4.0D+00*GRDAA
C
         CALL RKS_X_PW91(IDERIV,NPT,RHOA1,SIGMAAA1,ZK,VRHOA,
     *                   VSIGMAAA,V2RHOA2,V2RHOASIGMAAA,V2SIGMAAA2)
C
         WAA = VSIGMAAA
C
         VXCA1 = VXCA1 + VRHOA * CPW91
         DUMAX = DUMAX + 1.0D+00*(WAA*GRADXA) * CPW91
         DUMAY = DUMAY + 1.0D+00*(WAA*GRADYA) * CPW91
         DUMAZ = DUMAZ + 1.0D+00*(WAA*GRADZA) * CPW91
C
         XPW91 = XPW91 + ZK * FTOTWT * CPW91
      END IF
C
      IF(UROHF) THEN
C
         IDERIV   = 1
         NPT      = 1
         RHOA1    = ROA
         SIGMAAA1 = GRDAA
         RHOB1    = ROB
         SIGMABB1 = GRDBB
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
              WAA = VSIGMAAA
              WBB = VSIGMABB
C
            VXCA1 = VXCA1 + VRHOA * CPW91
            DUMAX = DUMAX + 2.0D+00*(WAA*GRADXA) * CPW91
            DUMAY = DUMAY + 2.0D+00*(WAA*GRADYA) * CPW91
            DUMAZ = DUMAZ + 2.0D+00*(WAA*GRADZA) * CPW91
C
            VXCB1 = VXCB1 + VRHOB * CPW91
            DUMBX = DUMBX + 2.0D+00*(WBB*GRADXB) * CPW91
            DUMBY = DUMBY + 2.0D+00*(WBB*GRADYB) * CPW91
            DUMBZ = DUMBZ + 2.0D+00*(WBB*GRADZB) * CPW91
C
            XPW91 = XPW91 + ZK * FTOTWT * CPW91
C
      END IF
C
      RETURN
      END
C*MODULE DFTXCA  *DECK PZ81CF
      SUBROUTINE PZ81CF(RHFGVB,UROHF,ROA,ROB,FTOTWT,ECF,VXCA1,VXCB1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
C
C        INTERFACED BY ROBERTO PEVERATI AND KIM BALDRIDGE, U.ZURICH
C
      IF(RHFGVB) THEN
         IDERIV   = 1
         NPT      = 1
         RHOA1    = 2.0D+00*ROA
C
         CALL RKS_C_PZ81(IDERIV,NPT,RHOA1,ZK,VRHOA,V2RHOA2)
C
         VXCA1 = VXCA1 + VRHOA
         ECF = ECF + ZK * FTOTWT
      END IF
C
      IF(UROHF) THEN
         IDERIV   = 1
         NPT      = 1
         RHOA1    = ROA
         RHOB1    = ROB
C
         CALL UKS_C_PZ81(IDERIV,NPT,RHOA1,RHOB1,ZK,VRHOA,VRHOB,
     *                   V2RHOA2,V2RHOB2,V2RHOAB)
C
         VXCA1 = VXCA1 + VRHOA
         VXCB1 = VXCB1 + VRHOB
         ECF = ECF + ZK * FTOTWT
      END IF
C
      RETURN
      END
C
C*MODULE DFTXCA  *DECK VWN1CF
      SUBROUTINE VWN1CF(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA1,VXCB1,EVWN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF, RHFGVB
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
C
C     THE FOLLOWING IS THE BASIC SETUP FOR INTEGRATING AN LDA TYPE
C     CORRELATION FUNCTIONAL FROM THE DFT REPOSITORY.
C
      IDERIV= 1
      NPT   = 1
C
      IF(RHFGVB)THEN
         RHOA1 = 2.0D+00*ROA
         CALL RKS_C_VWN1RPA(IDERIV,NPT,RHOA1,ZK,VRHOA,V2RHOA2)
         VXCA1 = VXCA1 +     VRHOA*CVWN
         EVWN  = EVWN  + ZK*FTOTWT*CVWN
      END IF
C
      IF(UROHF) THEN
         RHOA1 = ROA
         RHOB1 = ROB
         CALL UKS_C_VWN1RPA(IDERIV,NPT,RHOA1,RHOB1,ZK,
     &                      VRHOA,VRHOB,V2RHOA2,V2RHOB2,V2RHOAB)
         VXCA1 = VXCA1 +     VRHOA*CVWN
         VXCB1 = VXCB1 +     VRHOB*CVWN
         EVWN  = EVWN  + ZK*FTOTWT*CVWN
      END IF
C
      RETURN
      END
C*MODULE DFTXCA  *DECK VWN5CF
      SUBROUTINE VWN5CF(RHFGVB,UROHF,ROA,ROB,FTOTWT,VXCA,VXCB,EVWN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL UROHF, RHFGVB
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00,
     *           FOUR=4.0D+00, SIX=6.0D+00, EIGHT=8.0D+00,
     *           ANINE=9.0D+00)
      PARAMETER (THIRD=ONE/THREE, FRTHRD=FOUR/THREE, SIXTH=ONE/SIX)
C                FACR = (THREE/(FOUR*PI))**SIXTH
C                FACH = FOUR/((TWO**THIRD-ONE)*ANINE)
      PARAMETER (FACR = 0.78762331789974325053D+00)
      PARAMETER (FACH = 1.70992093416136561756D+00)
      PARAMETER (FACG = ANINE/EIGHT)
      PARAMETER (FACDG= THREE/TWO)
C
C        THE A,B,C,X0 PARAMETERS OF THE FIT ARE AS FOLLOWS.
C        P,F,A MEAN PARAMAGNETIC, FERROMAGNETIC, AND STIFFNESS TERMS.
C        NOTE UNITS OF THE -A- CONSTANTS ARE RYDBERGS, SO DIVIDE BY 2.
C
      PARAMETER (VWNAP= 0.0621814D+00/TWO, VWNBP= 3.7274400D+00,
     *           VWNCP=12.9352000D+00, VWNX0P=-0.1049800D+00)
      PARAMETER (VWNAF= VWNAP/TWO, VWNBF= 7.0604200D+00,
     *           VWNCF=18.0578000D+00, VWNX0F=-0.3250000D+00)
      PARAMETER (VWNAA=-0.03377372D+00/TWO, VWNBA= 1.1310700D+00,
     *           VWNCA=13.0045000D+00, VWNX0A=-0.0047584D+00)
C
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
C
C     VWN CORRELATION FUNCTIONAL - FORMULA V
C     S.H.VOSKO, L.WILK, M.NUSAIR  CAN.J.PHYS. 58, 1200-11(1980)
C     THE NECESSARY FORMULAE ARE CLEARLY SPELLED OUT IN THE APPENDIX OF
C     B.G.JOHNSON, P.M.W.GILL, J.A.POPLE  J.CHEM.PHYS. 98, 5612-26(1993)
C     THE PARAMETERS APPEARING ABOVE CAN BE FOUND IN THE MOLPRO MANUAL,
C     ALTHOUGH A FEW ARE VISIBLE IN THE VWN PAPER ITSELF.
C
C     ON ENTRY,
C     ROA,B   .... SPIN DENSITIES
C     FTOTWT  .... VOLUME OF THIS GRID POINT
C
C     ON EXIT,
C     EVWN    .... CORRELATION ENERGY OF VWN5
C     VXCA,B  .... FIRST DERIVATIVES OF THE CORRELATION FUNCTIONAL
C
      IF (RHFGVB) THEN
        FACVWN = CVWN
        RHOA =ROA
        RHOB =ROB
        RHOT =RHOA+RHOB
        RHO16=RHOT**SIXTH
C
        Q    = SQRT(FOUR*VWNCP-VWNBP*VWNBP)
        XX0  = VWNX0P*VWNX0P+VWNX0P*VWNBP+VWNCP
        COEF1= VWNAP
        COEF2=-VWNAP*VWNBP*VWNX0P/XX0
        COEF3= VWNAP*(TWO*VWNBP/Q)*((VWNCP-VWNX0P*VWNX0P)/XX0)
C
        X    =FACR/RHO16
        XX   =X*X+X*VWNBP+VWNCP
        T1   =LOG(X*X/XX)
        T2   =LOG((X-VWNX0P)*(X-VWNX0P)/XX)
        T3   =ATAN(Q/(TWO*X+VWNBP))
        EPSC =COEF1*T1+COEF2*T2+COEF3*T3
        EVWN =EVWN+FACVWN*EPSC*RHOT*FTOTWT
C
        DEPSC= VWNAP*(ONE/X-(X/XX)*(ONE+VWNBP/(X-VWNX0P)))
        DUMA = (EPSC-X*DEPSC/THREE)
        DUMB = ZERO
        VXCA = VXCA + FACVWN*DUMA
        VXCB = VXCB + DUMB
      END IF
C
      IF (UROHF) THEN
        FACVWN=CVWN
        RHOA  =ROA
        RHOB  =ROB
        RHOT  =RHOA+RHOB
        RHO16 =RHOT**SIXTH
C
        QP    = SQRT(FOUR*VWNCP-VWNBP*VWNBP)
        QF    = SQRT(FOUR*VWNCF-VWNBF*VWNBF)
        QA    = SQRT(FOUR*VWNCA-VWNBA*VWNBA)
        XX0P  = VWNX0P*VWNX0P+VWNX0P*VWNBP+VWNCP
        XX0F  = VWNX0F*VWNX0F+VWNX0F*VWNBF+VWNCF
        XX0A  = VWNX0A*VWNX0A+VWNX0A*VWNBA+VWNCA
C
C          THE THREE POTENTIALS
C
        COEF1P= VWNAP
        COEF1F= VWNAF
        COEF1A= VWNAA
        COEF2P=-VWNAP*VWNBP*VWNX0P/XX0P
        COEF2F=-VWNAF*VWNBF*VWNX0F/XX0F
        COEF2A=-VWNAA*VWNBA*VWNX0A/XX0A
        COEF3P= VWNAP*(TWO*VWNBP/QP)*((VWNCP-VWNX0P*VWNX0P)/XX0P)
        COEF3F= VWNAF*(TWO*VWNBF/QF)*((VWNCF-VWNX0F*VWNX0F)/XX0F)
        COEF3A= VWNAA*(TWO*VWNBA/QA)*((VWNCA-VWNX0A*VWNX0A)/XX0A)
C
        X     =FACR/RHO16
        XXP   =X*X+X*VWNBP+VWNCP
        XXF   =X*X+X*VWNBF+VWNCF
        XXA   =X*X+X*VWNBA+VWNCA
        T1P   =LOG(X*X/XXP)
        T1F   =LOG(X*X/XXF)
        T1A   =LOG(X*X/XXA)
        T2P   =LOG((X-VWNX0P)*(X-VWNX0P)/XXP)
        T2F   =LOG((X-VWNX0F)*(X-VWNX0F)/XXF)
        T2A   =LOG((X-VWNX0A)*(X-VWNX0A)/XXA)
        T3P   =ATAN(QP/(TWO*X+VWNBP))
        T3F   =ATAN(QF/(TWO*X+VWNBF))
        T3A   =ATAN(QA/(TWO*X+VWNBA))
C
        EPSCP =COEF1P*T1P+COEF2P*T2P+COEF3P*T3P
        EPSCF =COEF1F*T1F+COEF2F*T2F+COEF3F*T3F
        EPSCA =COEF1A*T1A+COEF2A*T2A+COEF3A*T3A
C
C          DERIVATIVE OF THE THREE POTENTIALS
C
        XBP    = TWO*X+VWNBP
        XBF    = TWO*X+VWNBF
        XBA    = TWO*X+VWNBA
        XBQ2P  = XBP*XBP + QP*QP
        XBQ2F  = XBF*XBF + QF*QF
        XBQ2A  = XBA*XBA + QA*QA
        BRACKP = TWO/(X-VWNX0P) -XBP/XXP -FOUR*(TWO*VWNX0P+VWNBP)/XBQ2P
        BRACKF = TWO/(X-VWNX0F) -XBF/XXF -FOUR*(TWO*VWNX0F+VWNBF)/XBQ2F
        BRACKA = TWO/(X-VWNX0A) -XBA/XXA -FOUR*(TWO*VWNX0A+VWNBA)/XBQ2A
        CURLYP = TWO/X - XBP/XXP - FOUR*VWNBP/XBQ2P
        CURLYF = TWO/X - XBF/XXF - FOUR*VWNBF/XBQ2F
        CURLYA = TWO/X - XBA/XXA - FOUR*VWNBA/XBQ2A
        DEPSCP = VWNAP*(CURLYP - VWNBP*VWNX0P*BRACKP/XX0P)
        DEPSCF = VWNAF*(CURLYF - VWNBF*VWNX0F*BRACKF/XX0F)
        DEPSCA = VWNAA*(CURLYA - VWNBA*VWNX0A*BRACKA/XX0A)
C
C          SOME AUXILIARY COMBINATIONS
C
        HX  = FACH*(EPSCF-EPSCP)/EPSCA - ONE
        DHX = FACH*(DEPSCF-DEPSCP-(EPSCF-EPSCP)*DEPSCA/EPSCA)/EPSCA
C
        ZETA  = (RHOA-RHOB)/RHOT
        ZETA2 = ZETA *ZETA
        ZETA3 = ZETA2*ZETA
        ZETA4 = ZETA2*ZETA2
        GZETA = ((ONE+ZETA)**FRTHRD+(ONE-ZETA)**FRTHRD-TWO)*FACG
        DGZETA= ((ONE+ZETA)**THIRD -(ONE-ZETA)**THIRD     )*FACDG
C
C          TOTAL POTENTIAL, AND CONTRIBUTION TO ENERGY
C
        POT   = EPSCP + EPSCA*GZETA*(ONE+HX*ZETA4)
        EVWN  = EVWN + FACVWN*POT*RHOT*FTOTWT
C
C          CONTRIBUTION TO FUNCTIONAL DERIVATIVE
C
        DPOT1 = -(X/SIX)*(DEPSCP + DEPSCA*GZETA*(ONE+HX*ZETA4)
     *                           +  EPSCA*GZETA*DHX*ZETA4)
        DPOT2 = EPSCA*(DGZETA*(ONE+HX*ZETA4) + FOUR*GZETA*HX*ZETA3)
        DUMA  = POT + (DPOT1 + DPOT2*(ONE-ZETA))
        DUMB  = POT + (DPOT1 - DPOT2*(ONE+ZETA))
C
        VXCA  = VXCA + FACVWN*DUMA
        VXCB  = VXCB + FACVWN*DUMB
      ENDIF
C
      RETURN
      END
