C 25 MAR 00 - MWS - EIGENVALUE CUTOFF ADDED TO 11 FUNCTIONAL ROUTINES
C 10 MAR 99 - KRG - NEW MODULE FOR DFT FUNCTIONALS
C
C*MODULE DFTFUN  *DECK PWLOC
      SUBROUTINE PWLOC(R,RSA,ROPEN,RSB,S,UM,TM,WRK,IWRK,HISPIN,
     *                 RTX,RTY,RTZ,RX2,RY2,RZ2,NUM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL HISPIN
C
      DIMENSION R(NUM,NUM),RSA(NUM,NUM),
     *          RTX(NUM,NUM),RTY(NUM,NUM),RTZ(NUM,NUM),
     *          S(NUM,NUM),RSB(NUM,NUM),ROPEN(NUM,NUM),
     *          UM(NUM,NUM),TM(NUM,NUM),WRK(NUM,8),IWRK(NUM),
     *          RX2(NUM,NUM),RY2(NUM,NUM),RZ2(NUM,NUM)
C
      PARAMETER (PI=3.141592653589793238462643D+00)
      PARAMETER (ONESIX=0.166666666666666666666666666667D+00)
      PARAMETER (FORTHRD=1.33333333333333333333333333333D+00)
      PARAMETER (ONETHRD=0.33333333333333333333333333333D+00)
C
      L3=NUM*NUM
C
C  THIS IS THE LOCAL PW CORRELATION FUNCTIONAL
C  IT IS DESCRIBED IN PHYSICAL REVIEW B, 45, 13244 (1992),
C  JOHN P PERDEW AND YUE WANG.
C
      IF(.NOT. HISPIN) THEN
C
C     SPECIAL RHF CASE, SINCE MOST TERM SIMPLIFY A LOT
C     THE G(WIGGLE) TERM IS ZERO SINCE ALPHA-BETA=ZERO
C
        CALL DFTJAC(R,UM,NUM,TM,IWRK,WRK,S)
        EIGCUT=1.0D-09
C
        DO 830 I=1,NUM
        IF (R(I,1).LT.EIGCUT) R(I,1) = EIGCUT
        IF (R(I,1).LT.EIGCUT) THEN
          RSA(I,1)   = 0.0D+00
          R  (I,1)   = 0.0D+00
        ELSE
C         MULTIPLY BY TWO TO GET TOTAL DENSITY REPRESENTATION OF X
          DUM      = (3.0D+00/(8.0D+00*PI*R(I,1)))**ONESIX
          RSA(I,1) = PWLOC1(DUM)
          R  (I,1) = RSA(I,1)-(DUM*PWLOC4(DUM)*ONESIX)
        END IF
 830    CONTINUE
        CALL DFTTRN1(RSA,UM,TM,NUM)
        CALL DFTTRN1(R,UM,TM,NUM)
C
C       HIGH SPIN CASE
C
      ELSE
        DO 100 I=1,NUM
          DO 100 J=1,NUM
            RTX(J,I)=R(J,I)+ROPEN(J,I)
            RX2(J,I)=R(J,I)-ROPEN(J,I)
C
C           RTX IS THE TOTAL DENSITY AND RX2 IS THE DENSITY DIFFERENCE
C
 100    CONTINUE
        CALL DFTJAC(RTX,UM ,NUM,TM,IWRK,WRK,S)
        EIGCT2=1.0D-11
C
        DO 101 I=1,NUM
          IF(RTX(I,1) .LT. EIGCT2) THEN
             RY2(I,1)=0.0D+00
          ELSE
             RY2(I,1)=1.0D+00/RTX(I,1)
          END IF
 101    CONTINUE
        CALL DFTTRN1(RY2,UM,TM,NUM)
C
        DO 102 I=1,NUM
         DO 103 J=1,I
          RTY(J,I)=0.0D+00
          DO 104 M=1,NUM
           RTY(J,I)=RTY(J,I)+(RX2(M,J)*RY2(M,I)+RY2(M,J)*RX2(M,I))
 104      CONTINUE
          RTY(J,I)=RTY(J,I)*0.5D+00
          RTY(I,J)=RTY(J,I)
 103     CONTINUE
 102    CONTINUE
C
C       RTY IS THE DENSITY POLARIZATION
C
        CALL DFTJAC(RTY,RTZ,NUM,TM,IWRK,WRK,S)
        DO 105 I=1,NUM
          IF(RTY(I,1) .GE.  1.0D+00) RTY(I,1)=1.0D+00
          IF(RTY(I,1) .LE. -1.0D+00) RTY(I,1)=-1.0D+00
 105    CONTINUE
C       UM AND RTZ ARE TRANSPOSED EIGENVECTORS
        DO 110 I=1,NUM
          IF(RTX(I,1) .LE. 0.0D+00)THEN
             RSA(I,1)=0.0D+00
             RY2(I,1)=0.0D+00
             RX2(I,1)=0.0D+00
             RTX(I,1)=0.0D+00
          ELSE
C            RTX IS NOW X
             DUM     =(3.0D+00/(4.0D+00*PI*RTX(I,1)))**ONESIX
             RTX(I,1)=DUM
             RSA(I,1)=PWLOC1(DUM)
             RY2(I,1)=PWLOC3(DUM)
             RX2(I,1)=RY2(I,1)*PWLOC7(DUM)
          END IF
 110    CONTINUE
        CALL DFTTRN1(RSA,UM,TM,NUM)
        CALL DFTTRN1(RY2,UM,TM,NUM)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        DO 120 I=1,NUM
          RZ2(I,1)=(9.0D+00/8.0D+00)*((1.0D+00+RTY(I,1))**(FORTHRD)+
     *                  (1.0D+00-RTY(I,1))**(FORTHRD)-2.0D+00)
          S(I,1)=RZ2(I,1)*RTY(I,1)**4
 120    CONTINUE
        CALL DFTTRN1(RZ2,RTZ,TM,NUM)
        CALL DFTTRN1(S,RTZ,TM,NUM)
C
C  RSA IS EP(RHO), RY2 IS EA(RHO), RZ2 IS G(WIGGLE),
C  RX2 IS H(RHO)*EA(RHO), S IS G(WIGGLE)*WIGGLE**4
C
        DO 130 J=1,NUM
          DO 130 I=1,J
            DO 131 M=1,NUM
C       SYMETRIZE THE FINAL ENERGY CONTRIBUTION, SOME INDEXES REVERSED
            RSA(I,J)=RSA(I,J)+(RY2(M,I)*RZ2(M,J)+RY2(M,J)*RZ2(M,I)+
     *                         RX2(M,I)*S(M,J)  +RX2(M,J)*S(M,I))
     *                    *0.5D+00
 131        CONTINUE
            RSA(J,I)=RSA(I,J)
 130    CONTINUE
      CALL DCOPY(NUM*NUM,RSA,1,RSB,1)
C
C     NOW THE FOCK CONTRIBUTION FOR HISPIN
C
C     NOTE: RTX IS "X" AND RTY IS (A-B)/(A+B) AS EIGENVALUES
C     UM AND RTZ ARE THE CORRESPONDING TRANSPOSED EIGENVECTORS
C
        CALL DCOPY(NUM*NUM,RSA,1,R,1)
        DO 140 I=1,NUM
          IF(RTX(I,1) .LE. 0.0D+00)THEN
             WRK(I,1)=0.0D+00
             WRK(I,2)=0.0D+00
             WRK(I,3)=0.0D+00
          ELSE
             WRK(I,1)=PWLOC4(RTX(I,1))*RTX(I,1)
             WRK(I,2)=PWLOC6(RTX(I,1))*RTX(I,1)
             WRK(I,3)=PWLOC6(RTX(I,1))*PWLOC7(RTX(I,1))*RTX(I,1)+
     *                PWLOC3(RTX(I,1))*PWLOC8(RTX(I,1))*RTX(I,1)
          END IF
 140    CONTINUE
        DO 150 I=1,NUM
          WRK(I,4)=(9.0D+00/8.0D+00)*((1.0D+00+RTY(I,1))**(FORTHRD)+
     *                  (1.0D+00-RTY(I,1))**(FORTHRD)-2.0D+00)
          WRK(I,5)=WRK(I,4)*RTY(I,1)**4
 150    CONTINUE
C       RX2 IS THE X BASED PART RY2 IS THE WIGGLE BASED PART
        CALL DCOPY(NUM,WRK(1,1),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DAXPY(L3,-(1.0D+00/6.0D+00),RX2,1,R,1)
C
        CALL DCOPY(NUM,WRK(1,2),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,4),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 160 I=1,NUM
          DO 160 J=1,NUM
            DO 160 M=1,NUM
              R(I,J)=R(I,J)+(-1.0D+00/6.0D+00)*RX2(M,I)*RY2(M,J)
 160    CONTINUE
        CALL DCOPY(NUM,WRK(1,3),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,5),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 170 I=1,NUM
          DO 170 J=1,NUM
            DO 170 M=1,NUM
              R(I,J)=R(I,J)+(-1.0D+00/6.0D+00)*RX2(M,I)*RY2(M,J)
 170    CONTINUE
C
        DO 180 I=1,NUM
          IF(RTX(I,1) .LE. 0.0D+00)THEN
             WRK(I,1)=0.0D+00
             WRK(I,2)=0.0D+00
          ELSE
             WRK(I,1)=PWLOC3(RTX(I,1))
             WRK(I,2)=WRK(I,1)*PWLOC7(RTX(I,1))
          END IF
 180    CONTINUE
        DO 190 I=1,NUM
          IF(RTY(I,1) .EQ.  0.0D+00) THEN
          WRK(I,3)=0.0D+00
          WRK(I,4)=0.0D+00
          WRK(I,5)=0.0D+00
          WRK(I,6)=0.0D+00
          ELSE
          WRK(I,3)=(3.0D+00/2.0D+00)*RTY(I,1)*
     *        ((1.0D+00+RTY(I,1))**(ONETHRD)-
     *         (1.0D+00-RTY(I,1))**(ONETHRD))
          WRK(I,4)=WRK(I,3)*RTY(I,1)**4+4.0D+00*RTY(I,1)**4*
     *            ((9.0D+00/8.0D+00)*((1.0D+00+RTY(I,1))**(FORTHRD)+
     *               (1.0D+00-RTY(I,1))**(FORTHRD)-2.0D+00))
          WRK(I,5)=WRK(I,3)/RTY(I,1)
          WRK(I,6)=WRK(I,4)/RTY(I,1)
          END IF
  190   CONTINUE
        CALL DCOPY(NUM,WRK(1,1),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,3),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 200 I=1,NUM
          DO 200 J=1,NUM
            DO 200 M=1,NUM
              R(I,J)=R(I,J)-RX2(M,I)*RY2(M,J)
 200    CONTINUE
        CALL DCOPY(NUM,WRK(1,2),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,4),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 210 I=1,NUM
          DO 210 J=1,NUM
            DO 210 M=1,NUM
              R(I,J)=R(I,J)-RX2(M,I)*RY2(M,J)
 210    CONTINUE
C
C      NOW THE PART THAT IS -/+ DEPENDING ON ALPHA VS. BETA FOCK
C
        CALL DCOPY(NUM,WRK(1,1),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,5),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 220 I=1,NUM
          DO 220 J=1,NUM
            DUM=0.0D+00
            DO 221 M=1,NUM
              DUM=DUM+RX2(M,I)*RY2(M,J)
 221        CONTINUE
            ROPEN(I,J)=R(I,J)-DUM
            R    (I,J)=R(I,J)+DUM
 220    CONTINUE
        CALL DCOPY(NUM,WRK(1,2),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,6),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 230 I=1,NUM
          DO 230 J=1,NUM
            DUM=0.0D+00
            DO 231 M=1,NUM
              DUM=DUM+RX2(M,I)*RY2(M,J)
 231        CONTINUE
            ROPEN(I,J)=ROPEN(I,J)-DUM
            R    (I,J)=R(I,J)+DUM
 230    CONTINUE
C       SYMMETRIZE THE FOCK
        DO 240 I=1,NUM
          DO 240 J=1,I
            R(I,J)=(R(I,J)+R(J,I))/2.0D+00
            R(J,I)=R(I,J)
            ROPEN(I,J)=(ROPEN(I,J)+ROPEN(J,I))/2.0D+00
            ROPEN(J,I)=ROPEN(I,J)
 240    CONTINUE
      ENDIF
      END
C*MODULE DFTFUN  *DECK PWLOC1
      DOUBLE PRECISION FUNCTION PWLOC1(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE G(X) FUNCTION FROM PW'S LOCAL FUNCTIONAL
C     THIS IS E(X,0), SIMILAR TO VWN'S EP
C     THIS IS NOT THE RPA
C     NOTE THAT THIS IS X NOT RS
C     PARAMETER (P=1.0D+00)
C     PARAMETER ( A=0.031091D+00)
C     AA=A+A
C
      PARAMETER (AA=0.062182D+00)
      PARAMETER (A1=0.2137D+00)
      PARAMETER (B1=7.5957D+00)
      PARAMETER (B2=3.5876D+00)
      PARAMETER (B3=1.6382D+00)
      PARAMETER (B4=0.49294D+00)
C
      PWLOC1=-AA*(1.0D+00+A1*X**2)*LOG(1.0D+00 + 1.0D+00/
     * (AA*(B1*X+B2*X**2+B3*X**3+B4*X**4)))
      END
C*MODULE DFTFUN  *DECK PWLOC2
      DOUBLE PRECISION FUNCTION PWLOC2(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE G(X) FUNCTION FROM PW'S LOCAL FUNCTIONAL
C     THIS IS E(X,1), SIMILAR TO VWN'S EF
C     THIS IS NOT THE RPA
C     NOTE THAT THIS IS X NOT RS
C     PARAMETER (P=1.0D+00)
C     PARAMETER (A=0.015545D+00)
C     AA=A+A
C
      PARAMETER (AA=0.03109D+00)
      PARAMETER (A1=0.20548D+00)
      PARAMETER (B1=14.1189D+00)
      PARAMETER (B2=6.1977D+00)
      PARAMETER (B3=3.3662D+00)
      PARAMETER (B4=0.62517D+00)
      PWLOC2=-AA*(1.0D+00+A1*X**2)*LOG(1.0D+00 + 1.0D+00/
     * (AA*(B1*X+B2*X**2+B3*X**3+B4*X**4)))
      END
C*MODULE DFTFUN  *DECK PWLOC3
      DOUBLE PRECISION FUNCTION PWLOC3(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE G(X) FUNCTION FROM PW'S LOCAL FUNCTIONAL
C     THIS IS A(X), SIMILAR TO VWN'S EA
C     THIS IS NOT THE RPA
C     NOTE THAT THIS IS X NOT RS
C     PARAMETER (P=1.0D+00)
C     PARAMETER (A=0.016887D+00)
C     AA=A+A
C
      PARAMETER (AA=0.033774D+00)
      PARAMETER (A1=0.11125D+00)
      PARAMETER (B1=10.357D+00)
      PARAMETER (B2=3.6231D+00)
      PARAMETER (B3=0.88026D+00)
      PARAMETER (B4=0.49671D+00)
      PWLOC3=-(-AA*(1.0D+00+A1*X**2)*LOG(1.0D+00 + 1.0D+00/
     * (AA*(B1*X+B2*X**2+B3*X**3+B4*X**4))))
      END
C*MODULE DFTFUN  *DECK PWLOC4
      DOUBLE PRECISION FUNCTION PWLOC4(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE G'(X) FUNCTION FROM PW'S LOCAL FUNCTIONAL
C     THIS IS E(X,0), SIMILAR TO VWN'S EP
C     THIS IS NOT THE RPA
C     NOTE THAT THIS IS X NOT RS
C     PARAMETER (P =1.0D+00)
C     PARAMETER (A =0.031091D+00)
C     AA=A+A
C
      PARAMETER (AA=0.062182D+00)
      PARAMETER (A1=0.2137D+00)
      PARAMETER (B1=7.5957D+00)
      PARAMETER (B2=3.5876D+00)
      PARAMETER (B3=1.6382D+00)
      PARAMETER (B4=0.49294D+00)
      BIGX=B1*X + B2*X**2 + B3*X**3 + B4*X**4
      PWLOC4=2.0D+00*AA*A1*X*LOG(1.0D+00+1.0D+00/(AA*BIGX))-
     *       ((1.0D+00+A1*X**2)/(BIGX**2+BIGX))*
     *       (B1+2.0D+00*B2*X+3.0D+00*B3*X**2+4.0D+00*B4*X**3)
      END
C*MODULE DFTFUN  *DECK PWLOC5
      DOUBLE PRECISION FUNCTION PWLOC5(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE G'(X) FUNCTION FROM PW'S LOCAL FUNCTIONAL
C     THIS IS E(X,1), SIMILAR TO VWN'S EF
C     THIS IS NOT THE RPA
C     NOTE THAT THIS IS X NOT RS
C     PARAMETER (P =1.0D+00)
C     PARAMETER (A =0.015545D+00)
C     AA=A+A
C
      PARAMETER (AA=0.03109D+00)
      PARAMETER (A1=0.20548D+00)
      PARAMETER (B1=14.1189D+00)
      PARAMETER (B2=6.1977D+00)
      PARAMETER (B3=3.3662D+00)
      PARAMETER (B4=0.62517D+00)
      BIGX=B1*X + B2*X**2 + B3*X**3 + B4*X**4
      PWLOC5=2.0D+00*AA*A1*X*LOG(1.0D+00+1.0D+00/(AA*BIGX))-
     *       ((1.0D+00+A1*X**2)/(BIGX**2+BIGX))*
     *       (B1+2.0D+00*B2*X+3.0D+00*B3*X**2+4.0D+00*B4*X**3)
      END
C*MODULE DFTFUN  *DECK PWLOC6
      DOUBLE PRECISION FUNCTION PWLOC6(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE G'(X) FUNCTION FROM PW'S LOCAL FUNCTIONAL
C     THIS IS A(X), SIMILAR TO VWN'S EA
C     THIS IS NOT THE RPA
C     NOTE THAT THIS IS X NOT RS
C     PARAMETER ( P=1.0D+00)
C     PARAMETER ( A=0.016887D+00)
C     AA=A+A
C
      PARAMETER (AA=0.033774D+00)
      PARAMETER (A1=0.11125D+00)
      PARAMETER (B1=10.357D+00)
      PARAMETER (B2=3.6231D+00)
      PARAMETER (B3=0.88026D+00)
      PARAMETER (B4=0.49671D+00)
      BIGX=B1*X + B2*X**2 + B3*X**3 + B4*X**4
      PWLOC6=-(2.0D+00*AA*A1*X*LOG(1.0D+00+1.0D+00/(AA*BIGX))-
     *       ((1.0D+00+A1*X**2)/(BIGX**2+BIGX))*
     *       (B1+2.0D+00*B2*X+3.0D+00*B3*X**2+4.0D+00*B4*X**3))
      END
C*MODULE DFTFUN  *DECK PWLOC7
      DOUBLE PRECISION FUNCTION PWLOC7(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS H(X) FUNCTION FROM PWLOCAL
C
      PWLOC7=1.7099212624D+00*(PWLOC2(X)-PWLOC1(X))/PWLOC3(X)-1.0D+00
      END
C*MODULE DFTFUN  *DECK PWLOC8
      DOUBLE PRECISION FUNCTION PWLOC8(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE DERIVATIVE OF H(X) FUNCTION PWLOCAL
C
      A11=PWLOC1(X)
      A12=PWLOC2(X)
      A13=PWLOC3(X)
      A14=PWLOC4(X)
      A15=PWLOC5(X)
      A16=PWLOC6(X)
      PWLOC8=1.7099212624D+00*(A15-A14-A16*(A12-A11)/A13)/A13
      END
C*MODULE DFTFUN  *DECK VWN
      SUBROUTINE VWN(R,RSA,ROPEN,RSB,S,UM,TM,WRK,IWRK,HISPIN,
     *               RTX,RTY,RTZ,RX2,RY2,RZ2,NUM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL HISPIN
C
      DIMENSION R(NUM,NUM),RSA(NUM,NUM),
     *          RTX(NUM,NUM),RTY(NUM,NUM),RTZ(NUM,NUM),
     *          S(NUM,NUM),RSB(NUM,NUM),ROPEN(NUM,NUM),
     *          UM(NUM,NUM),TM(NUM,NUM),WRK(NUM,8),IWRK(NUM),
     *          RX2(NUM,NUM),RY2(NUM,NUM),RZ2(NUM,NUM)
C
      PARAMETER (PI=3.141592653589793238462643D+00)
      PARAMETER (ONESIX=0.166666666666666666666666666667D+00)
      PARAMETER (FORTHRD=1.33333333333333333333333333333D+00)
      PARAMETER (ONETHRD=0.33333333333333333333333333333D+00)
C
      L3=NUM*NUM
C
C  THIS IS A VWN CORRELATION FUNCTIONAL
C  THIS IS THE VERSION OF VWN THAT THE INVENTORS PREFER
C  IT IS CALLED VWN V (KEYWORD VWN5) IN G92DFT
C  IT IS DESCRIBED IN CAN. J. PHYS. 58, 1200, (1980),
C  S. H. VOSKO, L. WILK, C  AND M. NUSAIR.
C  SEE PAGE 1209 RIGHT BEFORE THE EQUATION WITH NO NUMBER.
C  POPLE, GILL AND JOHNSON FURTHER EXPLAIN IT IN
C  J CHEM PHYS 98, 5612 (1993).
C
C     SPECIAL RHF CASE, SINCE MOST TERM SIMPLIFY A LOT
C     THE G(WIGGLE) TERM IS ZERO SINCE ALPHA-BETA=ZERO
C
      IF(.NOT. HISPIN)THEN
        CALL DFTJAC(R,UM,NUM,TM,IWRK,WRK,S)
        EIGCUT=1.0D-09
        DO 830 I=1,NUM
        IF (R(I,1).LT.EIGCUT) R(I,1) = EIGCUT
        IF (R(I,1).LT.EIGCUT) THEN
          RSA(I,1)   = 0.0D+00
          R  (I,1)   = 0.0D+00
        ELSE
C         MULTIPLY BY TWO TO GET TOTAL DENSITY REPRESENTATION OF X
          DUM      = (3.0D+00/(8.0D+00*PI*R(I,1)))**(ONESIX)
          RSA(I,1) = VWNFUN1(DUM)
          R  (I,1) = RSA(I,1)-(DUM*VWNFUN4(DUM)/6.0D+00)
        END IF
 830    CONTINUE
        CALL DFTTRN1(RSA,UM,TM,NUM)
        CALL DFTTRN1(R,UM,TM,NUM)
C
C       HIGH SPIN CASE
C
      ELSE
        DO 100 I=1,NUM
          DO 100 J=1,NUM
            RTX(J,I)=R(J,I)+ROPEN(J,I)
            RX2(J,I)=R(J,I)-ROPEN(J,I)
C
C           RTX IS THE TOTAL DENSITY AND RX2 IS THE DENSITY DIFFERENCE
C
 100    CONTINUE
        CALL DFTJAC(RTX,UM ,NUM,TM,IWRK,WRK,S)
        EIGCT2=1.0D-11
C
        DO 101 I=1,NUM
          IF(RTX(I,1) .LT. EIGCT2) THEN
            RY2(I,1)=0.0D+00
          ELSE
            RY2(I,1)=1.0D+00/RTX(I,1)
          END IF
 101    CONTINUE
        CALL DFTTRN1(RY2,UM,TM,NUM)
C
        DO 102 I=1,NUM
         DO 103 J=1,I
          RTY(J,I)=0.0D+00
          DO 104 M=1,NUM
           RTY(J,I)=RTY(J,I)+(RX2(M,J)*RY2(M,I)+RY2(M,J)*RX2(M,I))
 104      CONTINUE
          RTY(J,I)=RTY(J,I)*0.5D+00
          RTY(I,J)=RTY(J,I)
 103     CONTINUE
 102    CONTINUE
C
C       RTY IS THE DENSITY POLARIZATION
C
        CALL DFTJAC(RTY,RTZ,NUM,TM,IWRK,WRK,S)
        DO 105 I=1,NUM
          IF(RTY(I,1) .GE.  1.0D+00) RTY(I,1)=1.0D+00
          IF(RTY(I,1) .LE. -1.0D+00) RTY(I,1)=-1.0D+00
 105    CONTINUE
C       UM AND RTZ ARE TRANSPOSED EIGENVECTORS
        DO 110 I=1,NUM
          IF(RTX(I,1) .LE. 0.0D+00)THEN
             RSA(I,1)=0.0D+00
             RY2(I,1)=0.0D+00
             RX2(I,1)=0.0D+00
             RTX(I,1)=0.0D+00
          ELSE
C            RTX IS NOW X
             DUM     =(3.0D+00/(4.0D+00*PI*RTX(I,1)))**(ONESIX)
             RTX(I,1)=DUM
             RSA(I,1)=VWNFUN1(DUM)
             RY2(I,1)=VWNFUN3(DUM)
             RX2(I,1)=RY2(I,1)*VWNFUN7(DUM)
          END IF
 110    CONTINUE
        CALL DFTTRN1(RSA,UM,TM,NUM)
        CALL DFTTRN1(RY2,UM,TM,NUM)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        DO 120 I=1,NUM
          RZ2(I,1)=(9.0D+00/8.0D+00)*((1.0D+00+RTY(I,1))**(FORTHRD)+
     *                (1.0D+00-RTY(I,1))**(FORTHRD)-2.0D+00)
          S(I,1)=RZ2(I,1)*RTY(I,1)**4
 120    CONTINUE
        CALL DFTTRN1(RZ2,RTZ,TM,NUM)
        CALL DFTTRN1(S,RTZ,TM,NUM)
C
C  RSA IS EP(RHO), RY2 IS EA(RHO), RZ2 IS G(WIGGLE),
C  RX2 IS H(RHO)*EA(RHO), S IS G(WIGGLE)*WIGGLE**4
C
        DO 130 J=1,NUM
          DO 130 I=1,J
            DO 131 M=1,NUM
C        SYMETRIZE THE FINAL ENERGY CONTRIBUTION, SOME INDEXES REVERSED
            RSA(I,J)=RSA(I,J)+(RY2(M,I)*RZ2(M,J)+RY2(M,J)*RZ2(M,I)+
     *                         RX2(M,I)*S(M,J)+RX2(M,J)*S(M,I))*0.5D+00
 131        CONTINUE
            RSA(J,I)=RSA(I,J)
 130    CONTINUE
        CALL DCOPY(NUM*NUM,RSA,1,RSB,1)
C
C     NOW THE FOCK CONTRIBUTION FOR HISPIN
C
C     NOTE: RTX IS "X" AND RTY IS (A-B)/(A+B) AS EIGENVALUES
C     UM AND RTZ ARE THE CORRESPONDING TRANSPOSED EIGENVECTORS
C
        CALL DCOPY(NUM*NUM,RSA,1,R,1)
        DO 140 I=1,NUM
          IF(RTX(I,1) .LE. 0.0D+00)THEN
             WRK(I,1)=0.0D+00
             WRK(I,2)=0.0D+00
             WRK(I,3)=0.0D+00
          ELSE
             WRK(I,1)=VWNFUN4(RTX(I,1))*RTX(I,1)
             WRK(I,2)=VWNFUN6(RTX(I,1))*RTX(I,1)
             WRK(I,3)=VWNFUN6(RTX(I,1))*VWNFUN7(RTX(I,1))*RTX(I,1)+
     *                VWNFUN3(RTX(I,1))*VWNFUN8(RTX(I,1))*RTX(I,1)
          END IF
 140    CONTINUE
        DO 150 I=1,NUM
          WRK(I,4)=(9.0D+00/8.0D+00)*((1.0D+00+RTY(I,1))**(FORTHRD)+
     *                  (1.0D+00-RTY(I,1))**(FORTHRD)-2.0D+00)
          WRK(I,5)=WRK(I,4)*RTY(I,1)**4
 150    CONTINUE
C       RX2 IS THE X BASED PART RY2 IS THE WIGGLE BASED PART
        CALL DCOPY(NUM,WRK(1,1),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DAXPY(L3,-(1.0D+00/6.0D+00),RX2,1,R,1)
C
        CALL DCOPY(NUM,WRK(1,2),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,4),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 160 I=1,NUM
          DO 160 J=1,NUM
            DO 160 M=1,NUM
              R(I,J)=R(I,J)+(-1.0D+00/6.0D+00)*RX2(M,I)*RY2(M,J)
 160    CONTINUE
        CALL DCOPY(NUM,WRK(1,3),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,5),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 170 I=1,NUM
          DO 170 J=1,NUM
            DO 170 M=1,NUM
              R(I,J)=R(I,J)+(-1.0D+00/6.0D+00)*RX2(M,I)*RY2(M,J)
 170    CONTINUE
C
        DO 180 I=1,NUM
          IF(RTX(I,1) .LE. 0.0D+00)THEN
             WRK(I,1)=0.0D+00
             WRK(I,2)=0.0D+00
          ELSE
             WRK(I,1)=VWNFUN3(RTX(I,1))
             WRK(I,2)=WRK(I,1)*VWNFUN7(RTX(I,1))
          END IF
 180    CONTINUE
        DO 190 I=1,NUM
          IF(RTY(I,1) .EQ.  0.0D+00) THEN
          WRK(I,3)=0.0D+00
          WRK(I,4)=0.0D+00
          WRK(I,5)=0.0D+00
          WRK(I,6)=0.0D+00
          ELSE
          WRK(I,3)=(3.0D+00/2.0D+00)*RTY(I,1)*
     *        ((1.0D+00+RTY(I,1))**(ONETHRD)-
     *         (1.0D+00-RTY(I,1))**(ONETHRD))
          WRK(I,4)=WRK(I,3)*RTY(I,1)**4+4.0D+00*RTY(I,1)**4*
     *            ((9.0D+00/8.0D+00)*((1.0D+00+RTY(I,1))**(FORTHRD)+
     *               (1.0D+00-RTY(I,1))**(FORTHRD)-2.0D+00))
          WRK(I,5)=WRK(I,3)/RTY(I,1)
          WRK(I,6)=WRK(I,4)/RTY(I,1)
          END IF
 190    CONTINUE
        CALL DCOPY(NUM,WRK(1,1),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,3),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 200 I=1,NUM
          DO 200 J=1,NUM
            DO 200 M=1,NUM
              R(I,J)=R(I,J)-RX2(M,I)*RY2(M,J)
 200    CONTINUE
        CALL DCOPY(NUM,WRK(1,2),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,4),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 210 I=1,NUM
          DO 210 J=1,NUM
            DO 210 M=1,NUM
              R(I,J)=R(I,J)-RX2(M,I)*RY2(M,J)
 210    CONTINUE
C
C      NOW THE PART THAT IS -/+ DEPENDING ON ALPHA VS. BETA FOCK
C
        CALL DCOPY(NUM,WRK(1,1),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,5),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 220 I=1,NUM
          DO 220 J=1,NUM
            DUM=0.0D+00
            DO 221 M=1,NUM
              DUM=DUM+RX2(M,I)*RY2(M,J)
 221        CONTINUE
            ROPEN(I,J)=R(I,J)-DUM
            R    (I,J)=R(I,J)+DUM
 220    CONTINUE
        CALL DCOPY(NUM,WRK(1,2),1,RX2,1)
        CALL DFTTRN1(RX2,UM,TM,NUM)
        CALL DCOPY(NUM,WRK(1,6),1,RY2,1)
        CALL DFTTRN1(RY2,RTZ,TM,NUM)
        DO 230 I=1,NUM
          DO 230 J=1,NUM
            DUM=0.0D+00
            DO 231 M=1,NUM
              DUM=DUM+RX2(M,I)*RY2(M,J)
 231        CONTINUE
            ROPEN(I,J)=ROPEN(I,J)-DUM
            R    (I,J)=R(I,J)+DUM
 230    CONTINUE
C       SYMMETRIZE THE FOCK
        DO 240 I=1,NUM
          DO 240 J=1,I
            R(I,J)=(R(I,J)+R(J,I))/2.0D+00
            R(J,I)=R(I,J)
            ROPEN(I,J)=(ROPEN(I,J)+ROPEN(J,I))/2.0D+00
            ROPEN(J,I)=ROPEN(I,J)
 240    CONTINUE
      ENDIF
      END
C*MODULE DFTFUN  *DECK VWNFUN1
      DOUBLE PRECISION FUNCTION VWNFUN1(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE EP(X) FUNCTION FROM VWN
C
      PARAMETER (A=0.0310907D+00)
      PARAMETER (X0=-0.10498D+00)
      PARAMETER (B=3.72744D+00)
      PARAMETER (C=12.9352D+00)
C     Q=SQRT(4.0D+00*C-B*B)
      PARAMETER (Q=6.15199081975907980D+00)
C     Q2=2/Q
      PARAMETER (Q2=0.325098014382004974D+00)
      THETAN=Q2*ATAN(1.0D+00/(Q2*X+B/Q))
      VWNFUN1=A*(LOG(X*X/(X*X+B*X+C))+B*THETAN
     *   -(B*X0/(X0*X0+B*X0+C))*(LOG((X-X0)*(X-X0)/(X*X+B*X+C))
     *   +(B+2.0D+00*X0)*THETAN))
      END
C*MODULE DFTFUN  *DECK VWNFUN2
      DOUBLE PRECISION FUNCTION VWNFUN2(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE EF(X) FUNCTION FROM VWN
C
      PARAMETER (A=0.01553535D+00)
      PARAMETER (X0=-0.325D+00)
      PARAMETER (B=7.06042D+00)
      PARAMETER (C=18.0578D+00)
C     Q=SQRT(4.0D+00*C-B*B)
      PARAMETER (Q=4.73092690956011364D+00)
C     Q2=2/Q
      PARAMETER (Q2=0.422750137178923791D+00)
      THETAN=Q2*ATAN(1.0D+00/(Q2*X+B/Q))
      VWNFUN2=A*(LOG(X*X/(X*X+B*X+C))+B*THETAN
     *   -(B*X0/(X0*X0+B*X0+C))*(LOG((X-X0)*(X-X0)/(X*X+B*X+C))
     *   +(B+2.0D+00*X0)*THETAN))
      END
C*MODULE DFTFUN  *DECK VWNFUN3
      DOUBLE PRECISION FUNCTION VWNFUN3(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE EA(X) FUNCTION FROM VWN
C     ALSO CALLED A(X)
C
      PARAMETER (A=-0.01688686394D+00)
      PARAMETER (X0=-0.00475840D+00)
      PARAMETER (B=1.13107D+00)
      PARAMETER (C=13.0045D+00)
C     Q=SQRT(4.0D+00*C-B*B)
C     THETAN=(2.0D+00/Q)*ATAN(Q/(2.0D+00*X+B))
C     PARAMETER (Q=7.12310891781811772D+00)
C     Q2=2/Q
      PARAMETER (Q2=0.280776276633520949D+00)
C     BQ=B/Q
      PARAMETER (BQ=0.158788811605938274D+00)
      THETAN=Q2*ATAN(1.0D+00/(Q2*X+BQ))
      VWNFUN3=A*(LOG(X*X/(X*X+B*X+C))+B*THETAN
     *   -(B*X0/(X0*X0+B*X0+C))*(LOG((X-X0)*(X-X0)/(X*X+B*X+C))
     *   +(B+2.0D+00*X0)*THETAN))
      END
C*MODULE DFTFUN  *DECK VWNFUN4
      DOUBLE PRECISION FUNCTION VWNFUN4(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE DERIVATIVE OF THE EP(X) FUNCTION FROM VWN
C
      PARAMETER (A=0.0310907D+00)
      PARAMETER (X0=-0.10498D+00)
      PARAMETER (B=3.72744D+00)
      PARAMETER (C=12.9352D+00)
C     Q2=(4.0D+00*C-B*B)
      PARAMETER (Q2=37.8469910463999994D+00)
      THEQ=(2.0D+00*X+B)**2+Q2
      VWNFUN4=A*(2.0D+00/X-((X+X+B)/(X*X+B*X+C))-(4.0D+00*B/THEQ)-
     *    (B*X0/((X0*X0+B*X0+C)))*(
     *     2.0D+00/(X-X0)-((X+X+B)/(X*X+B*X+C))-4.0D+00*(X0+X0+B)/THEQ))
      END
C*MODULE DFTFUN  *DECK VWNFUN5
      DOUBLE PRECISION FUNCTION VWNFUN5(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE DERIVATIVE OF THE EF(X) FUNCTION FROM VWN
C
      PARAMETER (A=0.01553535D+00)
      PARAMETER (X0=-0.325D+00)
      PARAMETER (B=7.06042D+00)
      PARAMETER (C=18.0578D+00)
C     Q2=(4.0D+00*C-B*B)
      PARAMETER (Q2=22.3816694236000089D+00)
      THEQ=(2.0D+00*X+B)**2+Q2
      VWNFUN5=A*(2.0D+00/X-((X+X+B)/(X*X+B*X+C))-(4.0D+00*B/THEQ)-
     *    (B*X0/((X0*X0+B*X0+C)))*(
     *     2.0D+00/(X-X0)-((X+X+B)/(X*X+B*X+C))-4.0D+00*(X0+X0+B)/THEQ))
      END
C*MODULE DFTFUN  *DECK VWNFUN6
      DOUBLE PRECISION FUNCTION VWNFUN6(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE DERIVATIVE OF THE EA(X) FUNCTION FROM VWN
C     ALSO CALLED A(X)
C
      PARAMETER (A=-0.01688686394D+00)
      PARAMETER (X0=-0.00475840D+00)
      PARAMETER (B=1.13107D+00)
      PARAMETER (C=13.0045D+00)
C     Q2=(4.0D+00*C-B*B)
      PARAMETER (Q2=50.7386806550999978D+00)
      THEQ=(2.0D+00*X+B)**2+Q2
      VWNFUN6=A*(2.0D+00/X-((X+X+B)/(X*X+B*X+C))-(4.0D+00*B/THEQ)-
     *    (B*X0/((X0*X0+B*X0+C)))*(
     *     2.0D+00/(X-X0)-((X+X+B)/(X*X+B*X+C))-4.0D+00*(X0+X0+B)/THEQ))
      END
C*MODULE DFTFUN  *DECK VWNFUN7
      DOUBLE PRECISION FUNCTION VWNFUN7(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS H(X) FUNCTION FROM VWN
C
      VWNFUN7=1.7099212624D+00*(VWNFUN2(X)-VWNFUN1(X))/VWNFUN3(X)
     *        - 1.0D+00
      END
C*MODULE DFTFUN  *DECK VWNFUN8
      DOUBLE PRECISION FUNCTION VWNFUN8(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS IS THE DERIVATIVE OF H(X) FUNCTION FROM VWN
C
      DF1=VWNFUN1(X)
      DF2=VWNFUN2(X)
      DF3=VWNFUN3(X)
      DF4=VWNFUN4(X)
      DF5=VWNFUN5(X)
      DF6=VWNFUN6(X)
      VWNFUN8=1.7099212624D+00*(DF5-DF4-DF6*(DF2-DF1)/DF3)/DF3
      END
C*MODULE DFTFUN  *DECK WIGEXP
      SUBROUTINE WIGEXP(R,RSA,ROPEN,RSB,S,UM,TM,WRK,IWRK,HISPIN,NUM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL HISPIN
C
      DIMENSION R(NUM,NUM),RSA(NUM,NUM),ROPEN(NUM,NUM),RSB(NUM,NUM),
     *          S(NUM,NUM),UM(NUM,NUM),TM(NUM,NUM),WRK(NUM,8),IWRK(NUM)
C
      PARAMETER(A0=0.906D+00)
      PARAMETER(TINYK=2.1987D-02)
      PARAMETER(ONETHRD=0.33333333333333333333333333D+00)
C
      L3=NUM*NUM
C
C  THIS FUNCTIONAL HAS NO DEPENDENCE ON SPIN POLARIZATION (A-B)
C  THEREFORE THE ENERGY COMPONENTS ARE DONE TOGETHER
C  THIS IS A LOCAL EXHANGE-CORRELATION FUNCTIONAL "WIGNER"
C  PHYSICAL REVIEW A, VOL 46, NUM 9, R5320,
C  QUINGSHENG ZHAO AND ROBERT G. PARR
C
C  THIS IS A MODIFIED VERSION TO AVOID A SINGULARITY AS DESCRIBED IN
C  J. CHEM. PHYS., VOL 99, NUM 11, 8765 (1993), GREGORY J. LAMING,
C  VOLKER TERMATH, NICHOLAS C. HANDY
C
      EIGCUT=1.0D-09
      IF(HISPIN)THEN
        CALL VADD (R,1,ROPEN,1,RSA,1,L3)
C
C       RSA IS NOW TOTAL DENSITY
C
C       RSA AND RSB AT THE SAME TIME-----------------------
C
        CALL DFTJAC(RSA,UM,NUM,TM,IWRK,WRK,S)
C
        DO 930 I=1,NUM
        IF (RSA(I,1).LT.EIGCUT) RSA(I,1) = EIGCUT
        IF (RSA(I,1).LT.EIGCUT) THEN
          RSA(I,1)   = 0.0D+00
            R(I,1)   = 0.0D+00
        ELSE
          DUM        = RSA(I,1)**ONETHRD
            R(I,1) = -(A0/3.0D+00)*EXP(-TINYK*DUM)
     *               *(-TINYK*DUM*DUM+4.0D+00*DUM)
          RSA(I,1)=-A0*(DUM*EXP(-TINYK*DUM))
        END IF
 930    CONTINUE
        CALL DFTTRN1(RSA,UM,TM,NUM)
        CALL DCOPY(NUM*NUM,RSA,1,RSB,1)
        CALL DFTTRN1(R,UM,TM,NUM)
        CALL DCOPY(NUM*NUM,R,1,ROPEN,1)
C
C
      ELSE
        CALL DFTJAC(R,UM,NUM,TM,IWRK,WRK,S)
        DO 830 I=1,NUM
        IF (R(I,1).LT.EIGCUT) R(I,1) = EIGCUT
        IF (R(I,1).LT.EIGCUT) THEN
          RSA(I,1)   = 0.0D+00
          R  (I,1)   = 0.0D+00
        ELSE
C         MULTIPLY BY TWO TO GET TOTAL DENSITY
          DUM      = (2.0D+00*R(I,1))**ONETHRD
          RSA(I,1) =-A0*(DUM*EXP(-TINYK*DUM))
          R(I,1)   = -(A0/3.0D+00)*EXP(-TINYK*DUM)
     *               *(-TINYK*DUM*DUM+4.0D+00*DUM)
        END IF
 830    CONTINUE
        CALL DFTTRN1(R  ,UM,TM,NUM)
        CALL DFTTRN1(RSA,UM,TM,NUM)
      END IF
      END
C*MODULE DFTFUN  *DECK WSDFT
      SUBROUTINE WSDFT(R,RSA,ROPEN,RSB,S,UM,TM,WRK,IWRK,HISPIN,NUM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL HISPIN
C
      DIMENSION R(NUM,NUM),RSA(NUM,NUM),ROPEN(NUM,NUM),RSB(NUM,NUM),
     *          S(NUM,NUM),UM(NUM,NUM),TM(NUM,NUM),WRK(NUM,8),IWRK(NUM)
C
      PARAMETER(A0=0.93222D+00)
      PARAMETER(TINYK=9.47362D-03)
      PARAMETER(ONETHRD=0.33333333333333333333333333333D+00)
C
      L3=NUM*NUM
C
C  THIS FUNCTIONAL HAS NO DEPENDENCE ON SPIN POLARIZATION (A-B)
C  THEREFORE THE ENERGY COMPONENTS ARE DONE TOGETHER
C  THIS IS A LOCAL EXHANGE-CORRELATION FUNCTIONAL "WIGNER SCALED"
C  PHYSICAL REVIEW A, VOL 46, NUM 9, R5320,
C  QUINGSHENG ZHAO AND ROBERT G. PARR
C
      EIGCUT=1.0D-09
      IF(HISPIN)THEN
        CALL VADD (R,1,ROPEN,1,RSA,1,L3)
C
C       RSA IS NOW TOTAL DENSITY !!!!!!
C
C       RSA AND RSB AT THE SAME TIME
C       R AND ROPEN AT THE SAME TIME (NOT SPIN POLARIZED, SO THEY
C            ARE THE SAME)
C
        CALL DFTJAC(RSA,UM,NUM,TM,IWRK,WRK,S)
        DO 930 I=1,NUM
        IF (RSA(I,1).LT.EIGCUT) RSA(I,1) = EIGCUT
        IF (RSA(I,1).LT.EIGCUT) THEN
          RSA(I,1)   = 0.0D+00
            R(I,1)   = 0.0D+00
        ELSE
               DUM = RSA(I,1)**ONETHRD
            R(I,1) = -(A0/3.0D+00)*DUM*(
     *   4.0D+00-5.0D+00*TINYK*DUM*LOG(1.0D+00+(1.0D+00/(TINYK*DUM)))
     *         +TINYK*DUM/(1.0D+00+TINYK*DUM))
          RSA(I,1)=-A0*DUM*(1-TINYK*DUM*LOG(1+(1.0D+00/(TINYK*DUM))))
        END IF
 930    CONTINUE
        CALL DFTTRN1(RSA,UM,TM,NUM)
        CALL DCOPY(NUM*NUM,RSA,1,RSB,1)
        CALL DFTTRN1(R,UM,TM,NUM)
        CALL DCOPY(NUM*NUM,R,1,ROPEN,1)
C
      ELSE
C
        CALL DFTJAC(R,UM,NUM,TM,IWRK,WRK,S)
        DO 830 I=1,NUM
        IF (R(I,1).LT.EIGCUT) R(I,1) = EIGCUT
        IF (R(I,1).LT.EIGCUT) THEN
          RSA(I,1)   = 0.0D+00
          R  (I,1)   = 0.0D+00
        ELSE
C         MULTIPLY BY TWO TO GET TOTAL DENSITY
          DUM = (2.0D+00*R(I,1))**ONETHRD
          RSA(I,1)=-A0*DUM*
     *               (1-TINYK*DUM*LOG(1.0D+00+(1.0D+00/(TINYK*DUM))))
          R(I,1) = -(A0/3.0D+00)*DUM*(
     *   4.0D+00-5.0D+00*TINYK*DUM*LOG(1.0D+00+(1.0D+00/(TINYK*DUM)))+
     *      TINYK*DUM/(1.0D+00+TINYK*DUM))
        END IF
 830    CONTINUE
        CALL DFTTRN1(R  ,UM,TM,NUM)
        CALL DFTTRN1(RSA,UM,TM,NUM)
      END IF
      END
C*MODULE DFTFUN  *DECK WIGNER
      SUBROUTINE WIGNER(R,RSA,ROPEN,RSB,S,UM,TM,WRK,IWRK,HISPIN,NUM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL HISPIN
C
      DIMENSION R(NUM,NUM),RSA(NUM,NUM),ROPEN(NUM,NUM),RSB(NUM,NUM),
     *          S(NUM,NUM),UM(NUM,NUM),TM(NUM,NUM),WRK(NUM,8),IWRK(NUM)
C
      PARAMETER(A0=0.906D+00)
      PARAMETER(TINYK=2.1987D-02)
      PARAMETER(ONETHRD=0.33333333333333333333333333D+00)
C
      L3=NUM*NUM
C
C  THIS FUNCTIONAL HAS NO DEPENDENCE ON SPIN POLARIZATION (A-B)
C  THEREFORE THE ENERGY COMPONENTS ARE DONE TOGETHER
C  THIS IS A LOCAL EXHANGE-CORRELATION FUNCTIONAL "WIGNER"
C  PHYSICAL REVIEW A, VOL 46, NUM 9, R5320,
C  QUINGSHENG ZHAO AND ROBERT G. PARR
C
      EIGCUT=1.0D-09
      IF(HISPIN) THEN
        CALL VADD (R,1,ROPEN,1,RSA,1,L3)
C
C       RSA IS NOW TOTAL DENSITY
C
C       RSA AND RSB AT THE SAME TIME-----------------------
C
        CALL DFTJAC(RSA,UM,NUM,TM,IWRK,WRK,S)
        DO 930 I=1,NUM
        IF (RSA(I,1).LT.EIGCUT) RSA(I,1) = EIGCUT
        IF (RSA(I,1).LT.EIGCUT) THEN
          RSA(I,1)   = 0.0D+00
          R  (I,1)   = 0.0D+00
        ELSE
          DUM = RSA(I,1)**ONETHRD
          R(I,1) = -(A0/3.0D+00)*(DUM/(1.0D+00+TINYK*DUM))*
     *              (4.0D+00-(TINYK*DUM/(1.0D+00+TINYK*DUM)))
          RSA(I,1)=-A0*(DUM/(1.0D+00+TINYK*DUM))
        END IF
 930    CONTINUE
        CALL DFTTRN1(RSA,UM,TM,NUM)
        CALL DCOPY(NUM*NUM,RSA,1,RSB,1)
        CALL DFTTRN1(R,UM,TM,NUM)
        CALL DCOPY(NUM*NUM,R,1,ROPEN,1)
      ELSE
        CALL DFTJAC(R,UM,NUM,TM,IWRK,WRK,S)
        DO 830 I=1,NUM
        IF (R(I,1).LE.1.0D+09) THEN
          RSA(I,1)   = 0.0D+00
          R  (I,1)   = 0.0D+00
        ELSE
C         MULTIPLY BY TWO TO GET TOTAL DENSITY
          DUM      = (2.0D+00*R(I,1))**ONETHRD
          RSA(I,1) =-A0*(DUM/(1.0D+00+TINYK*DUM))
          R(I,1)   = -(A0/3.0D+00)*(DUM/(1.0D+00+TINYK*DUM))*
     *              (4.0D+00-(TINYK*DUM/(1.0D+00+TINYK*DUM)))
        END IF
 830    CONTINUE
        CALL DFTTRN1(R  ,UM,TM,NUM)
        CALL DFTTRN1(RSA,UM,TM,NUM)
      END IF
      END
C*MODULE DFTFUN  *DECK XALPHA
      SUBROUTINE XALPHA(R,RSA,S,UM,TM,WRK,IWRK,NUM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION R(NUM,NUM),RSA(NUM,NUM),S(NUM,NUM),
     *          UM(NUM,NUM),TM(NUM,NUM),WRK(NUM,8),IWRK(NUM)
C
C         CONST=-((3/PI)**(1/3))*(2**(1/3))
      PARAMETER (CONST=-1.24070098179879995D+00)
      PARAMETER (ONETHRD=0.33333333333333333333333333D+00)
C
C     THIS IS THE TRADITIONAL LSDA,SLATER,X-ALPHA,LSD,ETC.
C     EXCHANGE FUNCTIONAL WITH ALPHA=2/3
C
      CALL DFTJAC(R,UM,NUM,TM,IWRK,WRK,S)
      EIGCUT=1.0D-09
      DO 830 I=1,NUM
      IF (R(I,1).LT.EIGCUT) R(I,1) = EIGCUT
      IF (R(I,1).LT.EIGCUT) THEN
        R(I,1)  = 0.0D+00
      ELSE
        R(I,1)  = CONST*R(I,1)**ONETHRD
      END IF
 830  CONTINUE
C
      CALL DFTTRN1(R,UM,TM,NUM)
      DO 800 I=1,NUM
        DO 800 J=1,NUM
          RSA(I,J)=0.75D+00*R(I,J)
 800  CONTINUE
C
C      R     MATRIX IS: <C  *(DENSITY)**1/3>
C      RSA   MATRIX IS: <COP*(DENSITY)**1/3>
C
      RETURN
      END
C*MODULE DFTFUN  *DECK CAMA
      SUBROUTINE CAMA(R,RSA,DX,DY,DZ,S,RX2,RY2,RZ2,RTX,RTY,RTZ,
     *                R3,UM,TM,WRK,IWRK,NUM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION R(NUM,NUM),RSA(NUM,NUM),
     *          DX(NUM,NUM),DY(NUM,NUM),DZ(NUM,NUM),
     *          S(NUM,NUM),RX2(NUM,NUM),RY2(NUM,NUM),RZ2(NUM,NUM),
     *          RTX(NUM,NUM),RTY(NUM,NUM),RTZ(NUM,NUM),
     *          R3(NUM,NUM),UM(NUM,NUM),TM(NUM,NUM),WRK(NUM,8),IWRK(NUM)
C
C        CONST=-((3./PI)**(1./3.))*(2.**(1./3.))
C        CONST2=0.5D+00/CONST
      PARAMETER (CONST=-1.24070098179879995D+00)
      PARAMETER (CONST2=-0.402997988504117388D+00)
      PARAMETER (B = 0.008D+00)
      PARAMETER (ONETHRD=0.3333333333333333333333333D+00)
C
      L3=NUM*NUM
C
C  THIS IS THE "POPULAR" BECKE88 GRADIENT CORRECTED EXCHANGE FUNCTIONAL
C  WITH "IMPROVED" VALUES FOR THE CONSTANTS
C  PHYSICAL REVIEW A, VOL 28, NUM 6, 3098 (1988), A. D. BECKE
C  J. CHEM. PHYS., VOL 99, NUM 11, 8765 (1993), GREGORY J. LAMING,
C  VOLKER TERMATH, NICHOLAS C. HANDY
C
C     START BY DOING LSDA EXCHANGE
C
      CALL DFTJAC(R,UM,NUM,TM,IWRK,WRK,S)
      EIGCUT=1.0D-09
      DO 830 I=1,NUM
      IF (R(I,1).LT.EIGCUT) R(I,1) = EIGCUT
      IF (R(I,1).LT.EIGCUT) THEN
        R3(I,1)  = 0.0D+00
        R(I,1)   = 0.0D+00
      ELSE
        DUM      = R(I,1)**ONETHRD
        R3(I,1)  = 1.0D+00/DUM
        R(I,1)   = CONST*DUM
      END IF
 830  CONTINUE
      CALL DFTTRN1(R,UM,TM,NUM)
      CALL DFTTRN1(R3,UM,TM,NUM)
      CALL DCOPY(L3,R,1,RSA,1)
C
C     GENERATE X
C
      DO 201 J=1,NUM
      DO 201 I=1,J
      RTX(I,J)=0.0D+00
      RTY(I,J)=0.0D+00
      RTZ(I,J)=0.0D+00
      DO 202 M=1,NUM
C
C      TAKE ADVANTAGE OF ANTI-SYMMETRIC DX AND SYMMETRIC R3 TO SPEED
C      UP MEMORY ACCESS
C       RTX(I,J)=RTX(I,J) -3.0D+00*(DX(I,M)*R3(M,J)-R3(I,M)*DX(M,J))
C
      RTX(I,J)=RTX(I,J) +(DX(M,I)*R3(M,J)+R3(M,I)*DX(M,J))
      RTY(I,J)=RTY(I,J) +(DY(M,I)*R3(M,J)+R3(M,I)*DY(M,J))
      RTZ(I,J)=RTZ(I,J) +(DZ(M,I)*R3(M,J)+R3(M,I)*DZ(M,J))
 202  CONTINUE
      RTX(I,J)=RTX(I,J)*3.0D+00
      RTY(I,J)=RTY(I,J)*3.0D+00
      RTZ(I,J)=RTZ(I,J)*3.0D+00
      RTX(J,I)=RTX(I,J)
      RTY(J,I)=RTY(I,J)
      RTZ(J,I)=RTZ(I,J)
 201  CONTINUE
C
C   RTZ=M[(D RHO/DZ)/(RHO**(4/3))]
C
      DO I=1,NUM
        DO J=1,I
          RX2(I,J)=0.0D+00
          DO M=1,NUM
            RX2(I,J)=RX2(I,J)+RTX(M,I)*RTX(M,J)
     *                       +RTY(M,I)*RTY(M,J)
     *                       +RTZ(M,I)*RTZ(M,J)
          END DO
          RX2(J,I)=RX2(I,J)
        END DO
      END DO
C
C  RX2 = X SQUARED
C
      DO 316 J=1,NUM
      DO 316 I=1,NUM
      RZ2(I,J)=0.0D+00
      DO 316 N=1,NUM
C      RZ2(I,J)=RZ2(I,J)+RTX(I,N)*DX(N,J)+RTY(I,N)*DY(N,J)+
C     *         RTZ(I,N)*DZ(N,J)
C      REVERSE RTX INDEX TO SPEED UP MEMORY ACCESS
      RZ2(I,J)=RZ2(I,J)+RTX(N,I)*DX(N,J)
     *                 +RTY(N,I)*DY(N,J)
     *                 +RTZ(N,I)*DZ(N,J)
C
C   RZ2(I,J)=[I*(D RHO/DZ)/(RHO**(4/3))*(DJ/DZ)] + X PART+ Y PART
C
 316  CONTINUE
      CALL DFTJAC(RX2,UM,NUM,TM,IWRK,WRK,S)
      EIGCT2=1.0D-11
C
      DO 300 I=1,NUM
      IF (RX2(I,1).LE.EIGCT2) THEN
         RX2(I,1) = 0.0D+00
         RTY(I,1) = 0.0D+00
         RTZ(I,1) = 0.0D+00
      ELSE
         RXX = SQRT(RX2(I,1))
         RXXXX=RX2(I,1)
C
C     RXX=X RX2(I,1)=X*X=RXXXX
C     NASTY TERM'S OF J.CHEM.PHYS.,98,5612(1993)
C
      RTY(I,1)=(12.0D+00*B**2*RXX*(RXX/SQRT(1.0D+00+
     *   RXXXX)-LOG(RXX+SQRT(RXXXX+1.0D+00)))
     *   -2.0D+00*B)/(1.0D+00+12.0D+00*B*RXX*LOG(RXX+
     *   SQRT(RXXXX+1.0D+00)))**2
      RTZ(I,1)=4.0D+00/3.0D+00*B*RXXXX*(1.0D+00-12.0D+00*B*
     *     RXXXX/SQRT(1.0D+00+RXXXX))/
     *     (1.0D+00+12.0D+00*B*RXX*LOG(RXX+SQRT(RXXXX+
     *     1.0D+00)))**2
      END IF
C
C   RTY=G'(X)/X    RTZ=(4/3)*[G(X)-X*G'(X)]   (CORRECTION PART ONLY)
C
 300  CONTINUE
      CALL DFTTRN1(RTY,UM,TM,NUM)
      CALL DFTTRN1(RTZ,UM,TM,NUM)
      DO 315 J=1,NUM
      DO 315 I=1,J
      RY2(I,J)=0.0D+00
      DO 314 M=1,NUM
C  ADD UP THE BECKE CORRECTION TO THE LSDA FOCK
C  FIRST PART IS M[(4/3)RHO**(1/3)*(G(X)-XG'(X))] WHICH IS THE PARTIAL
C  OF B WITH RESPECT TO RHO, IGNORING LSDA TERMS.
C  NOTE WELL:  [EQ A6 IS WRONG IN JCP]
C  SECOND PART IS 2*(I G'(X)/X M)(M*(D RHO/DZ)/(RHO**(4/3))*(DJ/DZ))
C  WHICH IS TWO TIMES THE PARTIAL OF B WITH RESPECT TO GAMMA TIMES
C  THE DERIVATIVE OF THE DENSITY TIMES THE DERIVATIVE OF I*J (THIS
C  LEADS TO TWO TERMS)  A X ABOVE FROM RZ2 AND A X BELOW FROM RTY
C  CANCEL (SEE EQ. 7 OF JCP,98,5612)
C
C     INDEXES ON R AND RTZ REVERSED FOR FASTER MEMORY ACCESS
      RY2(I,J)=RY2(I,J)+(CONST2*(R(M,I)*RTZ(M,J)+R(M,J)*RTZ(M,I)))
     *                 +(RTY(M,I)*RZ2(M,J)+RTY(M,J)*RZ2(M,I))
 314  CONTINUE
      RY2(J,I)=RY2(I,J)
 315  CONTINUE
      DO 900 J=1,NUM
      DO 900 I=1,NUM
      R(I,J) = R(I,J)+RY2(I,J)
 900  CONTINUE
C
C     GENERATE RSA----------------------------------------------------
C
      DO 340 I=1,NUM
      IF (RX2(I,1).EQ.0.0D+00) THEN
C        WE KNOW THAT IT'S NOT NEGATIVE FROM ABOVE
      ELSE
      RX2(I,1) = B*RX2(I,1)/(1.0D+00+12.0D+00*B*SQRT(RX2(I,1))
     *         *LOG(SQRT(RX2(I,1))+SQRT(RX2(I,1)+1.0D+00)))
C
      END IF
 340  CONTINUE
      CALL DFTTRN1(RX2,UM,TM,NUM)
      DO 650 J=1,NUM
      DO 650 I=1,J
      RTY(I,J)=0.0D+00
      DO 655 M=1,NUM
C
C     TO SYMMETRIZE THE BECKE ENERGY CORRECTION TO LSDA:
C               RSA*RX2=0.5*(RSA*RX2+RX2*RSA)/CONST
C     INDEXES ON RSA REVERSED FOR FASTER MEMORY ACCESS
C
      RTY(I,J)=RTY(I,J)+(RSA(M,I)*RX2(M,J)+RSA(M,J)*RX2(M,I))
 655  CONTINUE
      RTY(J,I)=RTY(I,J)
 650  CONTINUE
      DO 902 J=1,NUM
      DO 902 I=1,NUM
      RSA(I,J) = 0.75D+00*RSA(I,J)-CONST2*RTY(I,J)
 902  CONTINUE
      END
C*MODULE DFTFUN  *DECK BECKE88
      SUBROUTINE BECKE88(R,RSA,DX,DY,DZ,S,
     *                   RX2,RY2,RZ2,RTX,RTY,RTZ,R3,UM,TM,WRK,IWRK,NUM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION R(NUM,NUM),RSA(NUM,NUM),
     *          DX(NUM,NUM),DY(NUM,NUM),DZ(NUM,NUM),
     *          S(NUM,NUM),RX2(NUM,NUM),RY2(NUM,NUM),RZ2(NUM,NUM),
     *          RTX(NUM,NUM),RTY(NUM,NUM),RTZ(NUM,NUM),
     *          R3(NUM,NUM),UM(NUM,NUM),TM(NUM,NUM),WRK(NUM,8),IWRK(NUM)
C
C        CONST=-((3/PI)**(1/3))*(2**(1/3))
C        CONST2=0.5D+00/CONST
C        B6=6*B, BBD=6*B*B, B2=2*B, B43=3/3*B, BB8=8*B*B
C
      PARAMETER (ONETHRD=0.33333333333333333333D+00)
      PARAMETER (CONST=-1.24070098179879995D+00)
      PARAMETER (CONST2=-0.402997988504117388D+00)
      PARAMETER (B  =0.0042D+00)
      PARAMETER (B6 =0.0252D+00)
      PARAMETER (BB6=0.00010584D+00)
      PARAMETER (B2 =0.0084D+00)
      PARAMETER (B43=0.0056D+00)
      PARAMETER (BB8=0.00014112D+00)
C
      L3=NUM*NUM
C
C  THIS IS THE "POPULAR" BECKE88 GRADIENT CORRECTED EXCHANGE FUNCTIONAL
C  PHYSICAL REVIEW A, VOL 28, NUM 6, 3098 (1988), A. D. BECKE
C
C     START BY DOING LSDA EXCHANGE
C
      CALL DFTJAC(R,UM,NUM,TM,IWRK,WRK,S)
      EIGCUT=1.0D-09
      DO 830 I=1,NUM
      IF (R(I,1).LT.EIGCUT) R(I,1) = EIGCUT
      IF (R(I,1).LT.EIGCUT) THEN
        R3(I,1)  = 0.0D+00
        R(I,1)   = 0.0D+00
      ELSE
        DUM      = R(I,1)**ONETHRD
        R(I,1)   = CONST*DUM
        R3(I,1)  = 1.0D+00/DUM
      END IF
 830  CONTINUE
C
      CALL DFTTRN1(R,UM,TM,NUM)
      CALL DFTTRN1(R3,UM,TM,NUM)
      CALL DCOPY(L3,R,1,RSA,1)
C
C     GENERATE X
C
      DO 201 J=1,NUM
      DO 201 I=1,J
      RTX(I,J)=0.0D+00
      RTY(I,J)=0.0D+00
      RTZ(I,J)=0.0D+00
      DO 202 M=1,NUM
C
C      TAKE ADVANTAGE OF ANTI-SYMMETRIC DX AND SYMMETRIC R3 TO SPEED
C      UP MEMORY ACCESS
C       RTX(I,J)=RTX(I,J) -3.0D+00*(DX(I,M)*R3(M,J)-R3(I,M)*DX(M,J))
C
      RTX(I,J)=RTX(I,J) +(DX(M,I)*R3(M,J)+R3(M,I)*DX(M,J))
      RTY(I,J)=RTY(I,J) +(DY(M,I)*R3(M,J)+R3(M,I)*DY(M,J))
      RTZ(I,J)=RTZ(I,J) +(DZ(M,I)*R3(M,J)+R3(M,I)*DZ(M,J))
 202  CONTINUE
      RTX(I,J)=RTX(I,J)*3.0D+00
      RTY(I,J)=RTY(I,J)*3.0D+00
      RTZ(I,J)=RTZ(I,J)*3.0D+00
      RTX(J,I)=RTX(I,J)
      RTY(J,I)=RTY(I,J)
      RTZ(J,I)=RTZ(I,J)
 201  CONTINUE
C
C   RTZ=M[(D RHO/DZ)/(RHO**(4/3))]
C
      DO I=1,NUM
        DO J=1,I
          RX2(I,J)=0.0D+00
          DO M=1,NUM
            RX2(I,J)=RX2(I,J)+RTX(M,I)*RTX(M,J)
     *                       +RTY(M,I)*RTY(M,J)
     *                       +RTZ(M,I)*RTZ(M,J)
          END DO
          RX2(J,I)=RX2(I,J)
        END DO
      END DO
C
C  RX2 = X SQUARED
C
      DO 316 J=1,NUM
      DO 316 I=1,NUM
      RZ2(I,J)=0.0D+00
      DO 317 N=1,NUM
C      RZ2(I,J)=RZ2(I,J)+RTX(I,N)*DX(N,J)+RTY(I,N)*DY(N,J)+
C     *         RTZ(I,N)*DZ(N,J)
C      REVERSE RTX INDEX TO SPEED UP MEMORY ACCESS
      RZ2(I,J)=RZ2(I,J)+RTX(N,I)*DX(N,J)
     *                 +RTY(N,I)*DY(N,J)
     *                 +RTZ(N,I)*DZ(N,J)
C
C   RZ2(I,J)=[I*(D RHO/DZ)/(RHO**(4/3))*(DJ/DZ)] + X PART+ Y PART
C
 317  CONTINUE
 316  CONTINUE
      CALL DFTJAC(RX2,UM,NUM,TM,IWRK,WRK,S)
      EIGCT2=1.0D-11
C
      DO 300 I=1,NUM
      IF (RX2(I,1).LE.EIGCT2) THEN
         RX2(I,1) = 0.0D+00
         RTY(I,1) = 0.0D+00
         RTZ(I,1) = 0.0D+00
      ELSE
         RXXXX=RX2(I,1)
         RXX  =SQRT(RXXXX)
         RXX1 =SQRT(RXXXX+1.0D+00)
         RXLOG=LOG(RXX+RXX1)
         RXBIG=1.0D+00/(1.0D+00+B6*RXX*RXLOG)
C
C     RXX=X RX2(I,1)=X*X=RXXXX
C     NASTY TERM'S OF J.CHEM.PHYS.,98,5612(1993)
C
         RTY(I,1)=(BB6*(RXXXX/RXX1-RXX*RXLOG)-B2)*RXBIG*RXBIG
         RTZ(I,1)=CONST2*RXXXX*(B43-BB8*RXXXX/RXX1)*RXBIG*RXBIG
         RX2(I,1)=B*RXXXX*RXBIG
      END IF
C
C RTY=G'(X)/X  RTZ=CONST2*(4/3)*[G(X)-X*G'(X)] (CORRECTION PART ONLY)
C
 300  CONTINUE
      CALL DFTTRN1(RTY,UM,TM,NUM)
      CALL DFTTRN1(RTZ,UM,TM,NUM)
      CALL DFTTRN1(RX2,UM,TM,NUM)
      DO 315 J=1,NUM
      DO 315 I=1,J
      RY2(I,J)=0.0D+00
      DO 314 M=1,NUM
C
C  ADD UP THE BECKE CORRECTION TO THE LSDA FOCK
C  FIRST PART IS M[(4/3)RHO**(1/3)*(G(X)-XG'(X))] WHICH IS THE PARTIAL
C  OF B WITH RESPECT TO RHO, IGNORING LSDA TERMS.
C  NOTE WELL:  [EQ A6 IS WRONG IN JCP]
C  SECOND PART IS 2*(I G'(X)/X M)(M*(D RHO/DZ)/(RHO**(4/3))*(DJ/DZ))
C  WHICH IS TWO TIMES THE PARTIAL OF B WITH RESPECT TO GAMMA TIMES
C  THE DERIVATIVE OF THE DENSITY TIMES THE DERIVATIVE OF I*J (THIS LEADS
C  TO TWO TERMS)  A X ABOVE FROM RZ2 AND A X BELOW FROM RTY DO NOT
C  CANCEL THE X ABOVE DOES NOT HAVE ABSOLUTE VALUES ON IT, BECAUSE IT
C  IS THE DOT PRODUCT OF THE DENSITY GRADIENT WITH THE GRADIENT OF THE
C  ORBITALS(SEE EQ. 7 OF JCP,98,5612)
C
C     INDEXES ON R AND RTZ REVERSED FOR FASTER MEMORY ACCESS
      RY2(I,J)=RY2(I,J)+R(M,I)*RTZ(M,J)+R(M,J)*RTZ(M,I)
     *                 +RTY(M,I)*RZ2(M,J)+RTY(M,J)*RZ2(M,I)
 314  CONTINUE
C     DO THIS BELOW
C     RY2(J,I)=RY2(I,J)
 315  CONTINUE
      DO 900 J=1,NUM
      DO 900 I=1,J
      R(I,J) = R(I,J)+RY2(I,J)
      R(J,I) = R(I,J)
 900  CONTINUE
C
C     GENERATE RSA----------------------------------------------------
C
      DO 650 J=1,NUM
      DO 650 I=1,J
      RTY(I,J)=0.0D+00
      DO 655 M=1,NUM
C
C     TO SYMMETRIZE THE BECKE ENERGY CORRECTION TO LSDA:
C               RSA*RX2=0.5*(RSA*RX2+RX2*RSA)/CONST
C     INDEXES ON RSA REVERSED FOR FASTER MEMORY ACCESS
C
      RTY(I,J)=RTY(I,J)+(RSA(M,I)*RX2(M,J)+RSA(M,J)*RX2(M,I))
 655  CONTINUE
C     DO THIS BELOW
C     RTY(J,I)=RTY(I,J)
 650  CONTINUE
      DO 902 J=1,NUM
      DO 902 I=1,J
      RSA(I,J) = 0.75D+00*RSA(I,J)-CONST2*RTY(I,J)
      RSA(J,I) = RSA(I,J)
 902  CONTINUE
      END
C*MODULE DFTFUN  *DECK CAMB
      SUBROUTINE CAMB(R,RSA,DX,DY,DZ,S,RX2,RY2,RZ2,RTX,RTY,RTZ,
     *                R3,UM,TM,WRK,IWRK,RSAVE,NUM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION R(NUM,NUM),RSA(NUM,NUM),
     *          DX(NUM,NUM),DY(NUM,NUM),DZ(NUM,NUM),RSAVE(NUM,NUM),
     *          S(NUM,NUM),RX2(NUM,NUM),RY2(NUM,NUM),RZ2(NUM,NUM),
     *          RTX(NUM,NUM),RTY(NUM,NUM),RTZ(NUM,NUM),
     *          R3(NUM,NUM),UM(NUM,NUM),TM(NUM,NUM),WRK(NUM,8),IWRK(NUM)
C
C     CONST=-((3./PI)**(1./3.))*(2.**(1./3.))
      PARAMETER (CONST=-1.24070098179879995D+00)
      PARAMETER (B = 0.008D+00)
      PARAMETER (P=0.02D+00)
      PARAMETER (ONETHRD=0.33333333333333333333333333333D+00)
C
C  THIS IS THE "POPULAR" BECKE88 GRADIENT CORRECTED EXCHANGE FUNCTIONAL
C  WITH "IMPROVED" DEFINITION OF "LOCAL" EXCHANGE AND NEW CONSTANTS
C  THIS FUNCTIONAL INCLUDES SOME LOCAL CORRELATION
C  PHYSICAL REVIEW A, VOL 28, NUM 6, 3098 (1988), A. D. BECKE
C  J. CHEM. PHYS., VOL 99, NUM 11, 8765 (1993), GREGORY J. LAMING,
C  VOLKER TERMATH, NICHOLAS C. HANDY
C
C     START BY DOING LSDA EXCHANGE
C
      CALL DFTJAC(R,UM,NUM,TM,IWRK,WRK,S)
      EIGCUT=1.0D-09
C
      DO 830 I=1,NUM
      IF (R(I,1).LT.EIGCUT) R(I,1) = EIGCUT
      IF (R(I,1).LT.EIGCUT) THEN
        R3(I,1)   = 0.0D+00
        R(I,1)    = 0.0D+00
        RSA(I,1)  = 0.0D+00
        RSAVE(I,1)= 0.0D+00
      ELSE
        DUM       = R(I,1)**ONETHRD
        R3(I,1)   = 1.0D+00/DUM
        RSAVE(I,1)= DUM
        RSA(I,1)  = 0.75D+00*CONST*DUM*EXP(P*DUM)
        R(I,1)    = (CONST/3.0D+00)*EXP(P*R(I,1))*
     *               (P*R(I,1)*R(I,1)+4.0D+00*R(I,1))
      END IF
 830  CONTINUE
      CALL DFTTRN1(R,UM,TM,NUM)
      CALL DFTTRN1(RSAVE,UM,TM,NUM)
      CALL DFTTRN1(R3,UM,TM,NUM)
      CALL DFTTRN1(RSA,UM,TM,NUM)
C
C     GENERATE X
C
      DO 201 J=1,NUM
      DO 201 I=1,J
      RTX(I,J)=0.0D+00
      RTY(I,J)=0.0D+00
      RTZ(I,J)=0.0D+00
      DO 202 M=1,NUM
C
C      TAKE ADVANTAGE OF ANTI-SYMMETRIC DX AND SYMMETRIC R3 TO SPEED
C      UP MEMORY ACCESS
C       RTX(I,J)=RTX(I,J) -3.0D+00*(DX(I,M)*R3(M,J)-R3(I,M)*DX(M,J))
C
      RTX(I,J)=RTX(I,J) +(DX(M,I)*R3(M,J)+R3(M,I)*DX(M,J))
      RTY(I,J)=RTY(I,J) +(DY(M,I)*R3(M,J)+R3(M,I)*DY(M,J))
      RTZ(I,J)=RTZ(I,J) +(DZ(M,I)*R3(M,J)+R3(M,I)*DZ(M,J))
 202  CONTINUE
      RTX(I,J)=RTX(I,J)*3.0D+00
      RTY(I,J)=RTY(I,J)*3.0D+00
      RTZ(I,J)=RTZ(I,J)*3.0D+00
      RTX(J,I)=RTX(I,J)
      RTY(J,I)=RTY(I,J)
      RTZ(J,I)=RTZ(I,J)
 201  CONTINUE
C
C   RTZ=M[(D RHO/DZ)/(RHO**(4/3))]
C
      DO I=1,NUM
        DO J=1,I
          RX2(I,J)=0.0D+00
          DO M=1,NUM
            RX2(I,J)=RX2(I,J)+RTX(M,I)*RTX(M,J)
     *                       +RTY(M,I)*RTY(M,J)
     *                       +RTZ(M,I)*RTZ(M,J)
          END DO
          RX2(J,I)=RX2(I,J)
        END DO
      END DO
C
C  RX2 = X SQUARED
C
      DO 316 J=1,NUM
      DO 316 I=1,NUM
      RZ2(I,J)=0.0D+00
      DO 316 N=1,NUM
C      RZ2(I,J)=RZ2(I,J)+RTX(I,N)*DX(N,J)+RTY(I,N)*DY(N,J)+
C     *         RTZ(I,N)*DZ(N,J)
C      REVERSE RTX INDEX TO SPEED UP MEMORY ACCESS
      RZ2(I,J)=RZ2(I,J)+RTX(N,I)*DX(N,J)
     *                 +RTY(N,I)*DY(N,J)
     *                 +RTZ(N,I)*DZ(N,J)
C
C   RZ2(I,J)=[I*(D RHO/DZ)/(RHO**(4/3))*(DJ/DZ)] + X PART+ Y PART
C
 316  CONTINUE
      CALL DFTJAC(RX2,UM,NUM,TM,IWRK,WRK,S)
      EIGCT2=1.0D-11
C
      DO 300 I=1,NUM
      IF (RX2(I,1).LE.EIGCT2) THEN
         RX2(I,1) = 0.0D+00
         RTY(I,1) = 0.0D+00
         RTZ(I,1) = 0.0D+00
      ELSE
         RXX = SQRT(RX2(I,1))
         RXXXX=RX2(I,1)
C
C     RXX=X RX2(I,1)=X*X=RXXXX
C     NASTY TERM'S OF J.CHEM.PHYS.,98,5612(1993)
C
      RTY(I,1)=(6.0D+00*B**2*RXX*(RXX/SQRT(1.0D+00+
     *   RXXXX)-LOG(RXX+SQRT(RXXXX+1.0D+00)))
     *   -2.0D+00*B)/(1.0D+00+6.0D+00*B*RXX*LOG(RXX+
     *   SQRT(RXXXX+1.0D+00)))**2
      RTZ(I,1)=4.0D+00/3.0D+00*B*RXXXX*(1.0D+00-6.0D+00*B*
     *     RXXXX/SQRT(1.0D+00+RXXXX))/
     *     (1.0D+00+6.0D+00*B*RXX*LOG(RXX+SQRT(RXXXX+
     *     1.0D+00)))**2
      END IF
C
C   RTY=G'(X)/X    RTZ=(4/3)*[G(X)-X*G'(X)]   (CORRECTION PART ONLY)
C
 300  CONTINUE
      CALL DFTTRN1(RTY,UM,TM,NUM)
      CALL DFTTRN1(RTZ,UM,TM,NUM)
      DO 315 J=1,NUM
      DO 315 I=1,J
      RY2(I,J)=0.0D+00
      DO 314 M=1,NUM
C  ADD UP THE BECKE CORRECTION TO THE LSDA FOCK
C  FIRST PART IS M[(4/3)RHO**(1/3)*(G(X)-XG'(X))] WHICH IS THE PARTIAL
C  OF B WITH RESPECT TO RHO, IGNORING LSDA TERMS.
C  NOTE WELL: EQ A6 IS WRONG IN JCP
C  SECOND PART IS 2*(I G'(X)/X M)(M*(D RHO/DZ)/(RHO**(4/3))*(DJ/DZ))
C  WHICH IS TWO TIMES THE PARTIAL OF B WITH RESPECT TO GAMMA TIMES
C  THE DERIVATIVE OF THE DENSITY TIMES THE DERIVATIVE OF I*J (THIS LEADS
C  TO TWO TERMS)  A X ABOVE FROM RZ2 AND A X BELOW FROM RTY CANCEL
C  (SEE EQ. 7 OF JCP,98,5612)
C
C     INDEXES ON RSAVE AND RTZ REVERSED FOR FASTER MEMORY ACCESS
      RY2(I,J)=RY2(I,J)+
     *      (0.5D+00*(RSAVE(M,I)*RTZ(M,J)+RSAVE(M,J)*RTZ(M,I)))+
     *      (RTY(M,I)*RZ2(M,J)+RTY(M,J)*RZ2(M,I))
 314  CONTINUE
      RY2(J,I)=RY2(I,J)
 315  CONTINUE
      DO 900 J=1,NUM
      DO 900 I=1,NUM
      R(I,J) = R(I,J)+RY2(I,J)
 900  CONTINUE
C
C     GENERATE RSA----------------------------------------------------
C
      DO 340 I=1,NUM
      IF (RX2(I,1).EQ.0.0D+00) THEN
C        WE KNOW THAT IT'S NOT NEGATIVE FROM ABOVE
      ELSE
      RX2(I,1) = B*RX2(I,1)/(1.0D+00+6.0D+00*B*SQRT(RX2(I,1))
     *         *LOG(SQRT(RX2(I,1))+SQRT(RX2(I,1)+1.0D+00)))
C
      END IF
 340  CONTINUE
      CALL DFTTRN1(RX2,UM,TM,NUM)
      DO 650 J=1,NUM
      DO 650 I=1,J
      RTY(I,J)=0.0D+00
      DO 655 M=1,NUM
C
C     TO SYMMETRIZE THE BECKE ENERGY CORRECTION TO LSDA:
C               RSAVE*RX2=0.5*(RSAVE*RX2+RX2*RSAVE)
C     INDEXES ON RSA REVERSED FOR FASTER MEMORY ACCESS
C
      RTY(I,J)=RTY(I,J)+(RSAVE(M,I)*RX2(M,J)+RSAVE(M,J)*RX2(M,I))
 655  CONTINUE
      RTY(J,I)=RTY(I,J)
 650  CONTINUE
      DO 902 J=1,NUM
      DO 902 I=1,NUM
      RSA(I,J) = RSA(I,J)-(0.5D+00)*RTY(I,J)
 902  CONTINUE
      END
C*MODULE DFTFUN  *DECK DEPRISTO
      SUBROUTINE DEPRISTO(R,RSA,DX,DY,DZ,S,RX2,RY2,RZ2,RTX,RTY,RTZ,
     *                    R3,UM,TM,WRK,IWRK,NUM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION R(NUM,NUM),RSA(NUM,NUM),
     *          DX(NUM,NUM),DY(NUM,NUM),DZ(NUM,NUM),
     *          S(NUM,NUM),RX2(NUM,NUM),RY2(NUM,NUM),RZ2(NUM,NUM),
     *          RTX(NUM,NUM),RTY(NUM,NUM),RTZ(NUM,NUM),
     *          R3(NUM,NUM),UM(NUM,NUM),TM(NUM,NUM),WRK(NUM,8),IWRK(NUM)
C
C     CONST=-((3./PI)**(1./3.))*(2.**(1./3.))
C     CONST2=0.5/CONST
C
      PARAMETER (CONST=-1.24070098179879995D+00)
      PARAMETER (CONST2=-0.402997988504117388D+00)
      PARAMETER (B = 0.00132327D+00)
      PARAMETER (A1 = 0.861504D+00)
      PARAMETER (B1 = 0.044286D+00)
      PARAMETER (ONETHRD=0.33333333333333333333333333D+00)
C
      L3=NUM*NUM
C
C  THIS IS A GRADIENT CORRECTED EXCHANGE FUNCTIONAL
C  J. CHEM. PHYS. 86, 1425 (1987),ANDREW E. DEPRISTO AND JOEL E. KRESS
C
C     START BY DOING LSDA EXCHANGE
C
      CALL DFTJAC(R,UM,NUM,TM,IWRK,WRK,S)
      EIGCUT=1.0D-09
C
      DO 830 I=1,NUM
      IF (R(I,1).LT.EIGCUT) R(I,1) = EIGCUT
      IF (R(I,1).LT.EIGCUT) THEN
        R3(I,1)  = 0.0D+00
        R(I,1)   = 0.0D+00
      ELSE
        DUM      = R(I,1)**ONETHRD
        R3(I,1)  = 1.0D+00/DUM
        R(I,1)   = CONST*DUM
      END IF
 830  CONTINUE
      CALL DFTTRN1(R,UM,TM,NUM)
      CALL DFTTRN1(R3,UM,TM,NUM)
      CALL DCOPY(L3,R,1,RSA,1)
C
C     ----- GENERATE R -----
C
      DO 201 J=1,NUM
      DO 201 I=1,J
      RTX(I,J)=0.0D+00
      RTY(I,J)=0.0D+00
      RTZ(I,J)=0.0D+00
      DO 202 M=1,NUM
C
C      TAKE ADVANTAGE OF ANTI-SYMMETRIC DX AND SYMMETRIC R TO SPEED
C      UP MEMORY ACCESS
C       RTX(I,J)=RTX(I,J) -3.0D+00*(DX(I,M)*R3(M,J)-R3(I,M)*DX(M,J))
      RTX(I,J)=RTX(I,J) +(DX(M,I)*R3(M,J)+R3(M,I)*DX(M,J))
      RTY(I,J)=RTY(I,J) +(DY(M,I)*R3(M,J)+R3(M,I)*DY(M,J))
      RTZ(I,J)=RTZ(I,J) +(DZ(M,I)*R3(M,J)+R3(M,I)*DZ(M,J))
 202  CONTINUE
      RTX(I,J)=RTX(I,J)*3.0D+00
      RTY(I,J)=RTY(I,J)*3.0D+00
      RTZ(I,J)=RTZ(I,J)*3.0D+00
      RTX(J,I)=RTX(I,J)
      RTY(J,I)=RTY(I,J)
      RTZ(J,I)=RTZ(I,J)
 201  CONTINUE
C
      DO I=1,NUM
        DO J=1,I
          RX2(I,J)=0.0D+00
          DO M=1,NUM
            RX2(I,J)=RX2(I,J)+RTX(M,I)*RTX(M,J)
     *                       +RTY(M,I)*RTY(M,J)
     *                       +RTZ(M,I)*RTZ(M,J)
          END DO
          RX2(J,I)=RX2(I,J)
        END DO
      END DO
C
C     RX2 IS X SQUARED
C
      DO 316 J=1,NUM
      DO 316 I=1,NUM
      RZ2(I,J)=0.0D+00
      DO 316 N=1,NUM
C      RZ2(I,J)=RZ2(I,J)+RTX(I,N)*DX(N,J)+RTY(I,N)*DY(N,J)+
C     *         RTZ(I,N)*DZ(N,J)
C      REVERSE RTX INDEX TO SPEED UP MEMORY
      RZ2(I,J)=RZ2(I,J)+RTX(N,I)*DX(N,J)+RTY(N,I)*DY(N,J)+
     *         RTZ(N,I)*DZ(N,J)
 316  CONTINUE
      CALL DFTJAC(RX2,UM,NUM,TM,IWRK,WRK,S)
      EIGCT2=1.0D-11
C
      DO 300 I=1,NUM
      IF (RX2(I,1).LE.EIGCT2) THEN
         RX2(I,1) = 0.0D+00
         RTY(I,1) = 0.0D+00
         RTZ(I,1) = 0.0D+00
      ELSE
      RXX = SQRT(RX2(I,1))
      RXXXX=RX2(I,1)
C
C     NASTY TERM'S LIKE IN (J.CHEM.PHYS.,98,5612(1993))
C
      RTY(I,1)=B*(2.0D+00 + 3.0D+00*A1*RXX+ A1*B1*RXX**3)
     *        /(1.0D+00+B1*RXXXX)**2
      RTZ(I,1)=-4.0D+00/3.0D+00*B*(RXXXX*(B1*RXXXX
     *    -1.0D+00-2.0D+00*A1*RXX))/(1.0D+00+B1*RXXXX)**2
      END IF
 300  CONTINUE
      CALL DFTTRN1(RTY,UM,TM,NUM)
      CALL DFTTRN1(RTZ,UM,TM,NUM)
      DO 315 J=1,NUM
      DO 315 I=1,J
      RY2(I,J)=0.0D+00
      DO 314 M=1,NUM
C
C     TO SYMMETRIZE :   RY2=0.5*(R*RTZ+RTZ*R)/CONST
C      REVERSE INDEXES ON R FOR FASTER MEMORY ACCESS
C
      RY2(I,J)=RY2(I,J)+CONST2*(R(M,I)*RTZ(M,J)+R(M,J)*RTZ(M,I))
 314  CONTINUE
      RY2(J,I)=RY2(I,J)
 315  CONTINUE
      DO 318 J=1,NUM
      DO 318 I=1,J
      RTX(I,J)=0.0D+00
      DO 319 M=1,NUM
C     REVERSE INDEXES ON RTY FOR FASTER MEMORY ACCESS
      RTX(I,J)=RTX(I,J)+(RTY(M,I)*RZ2(M,J)+RTY(M,J)*RZ2(M,I))
 319  CONTINUE
      RTX(I,J)=RTX(I,J)*(-4.0D+00/3.0D+00)
      RTX(J,I)=RTX(I,J)
 318  CONTINUE
      DO 901 I=1,NUM
      DO 901 J=1,NUM
 901  R(I,J) = R(I,J)+RTX(I,J)+RY2(I,J)
C
C     GENERATE RSA----------------------------------------------------
C
      DO 340 I=1,NUM
      IF (RX2(I,1).LE.0.0D+00) THEN
        RX2(I,1) = 0.0D+00
      ELSE
C
C     EQ.(6) IN DEPRISTO'S PAPER
C
        RX2(I,1) = B*RX2(I,1)*(1.0D+00+A1*SQRT(RX2(I,1)))
     *         /(1.0D+00+B1*RX2(I,1))
C
      END IF
 340  CONTINUE
      CALL DFTTRN1(RX2,UM,TM,NUM)
C
      DO 650 J=1,NUM
      DO 650 I=1,J
      RTY(I,J)=0.0D+00
      DO 655 M=1,NUM
C     INDEXES ON RSA REVERSED FOR FASTER MEMORY ACCESS
      RTY(I,J)=RTY(I,J)+(RSA(M,I)*RX2(M,J)+RSA(M,J)*RX2(M,I))
 655  CONTINUE
      RTY(J,I)=RTY(I,J)
 650  CONTINUE
      DO 900 J=1,NUM
      DO 900 I=1,NUM
      RSA(I,J) = 0.75D+00*RSA(I,J)-CONST2*RTY(I,J)
 900  CONTINUE
      END
C*MODULE DFTFUN  *DECK LYP
      SUBROUTINE LYP(R,RSA,ROPEN,RSB,DX,DY,DZ,S,RX2,RY2,RZ2,
     *               RTX,RTY,RTZ,R3,UM,TM,WRK,IWRK,RSAVE,HISPIN,NUM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL HISPIN
C
      DIMENSION R(NUM,NUM),RSA(NUM,NUM),
     *          DX(NUM,NUM),DY(NUM,NUM),DZ(NUM,NUM),RSAVE(NUM,NUM),
     *          S(NUM,NUM),RSB(NUM,NUM),ROPEN(NUM,NUM),
     *          RX2(NUM,NUM),RY2(NUM,NUM),RZ2(NUM,NUM),
     *          RTX(NUM,NUM),RTY(NUM,NUM),RTZ(NUM,NUM),
     *          R3(NUM,NUM),UM(NUM,NUM),TM(NUM,NUM),WRK(NUM,8),IWRK(NUM)
C
C     AND NOW SOME PARAMETERS (LIKE ALL GOOD DFT FUNCTIONALS)  :-)
      PARAMETER (A=0.04918D+00)
      PARAMETER (B=0.132D+00)
C     PARAMETER (C=0.2533D+00)
      PARAMETER (D=0.349D+00)
C     PARAM1=-C/11
      PARAMETER (PARAM1=-0.02302727272727272727272727272727272D+00)
C     PARAM2=C+D
      PARAMETER (PARAM2=0.6023D+00)
C     PARAM3=C*D
      PARAMETER (PARAM3=0.0884017D+00)
C     PARAM4=-C/15
      PARAMETER (PARAM4=-0.016886666666666666666666666666667D+00)
C     PARAM5=(C-10D)/3
      PARAMETER (PARAM5=-1.0789D+00)
C     PARAM6=CD/3
      PARAMETER (PARAM6=0.02946723333333333333333333333333D+00)
C     PARAM7=-11/3
      PARAMETER (PARAM7=-3.666666666666666666666666666666667D+00)
C     PARAM8=-2CD/3
      PARAMETER (PARAM8=-0.05893446666666666666666666667D+00)
C     PARAM9=-(C+D)/3
      PARAMETER (PARAM9=-0.2007666666666666666666666667D+00)
C     PARAM10=-CDD/3
      PARAMETER (PARAM10=-0.010284064433D+00)
C     CF=(3.0D+00/10.0D+00)*(3.0D+00*PI*PI)**(2.0D+00/3.0D+00)
      PARAMETER (CF=2.87123400018819108D+00)
C     ABCF=A*B*CF*(2.0D+00**(11.0D+00/3.0D+00))
      PARAMETER (ABCF=0.236705143194386014D+00)
C
      L3=NUM*NUM
C
C  THIS IS THE LYP GRADIENT CORRECTED CORRELATION FUNCTIONAL
C  IT IS BASED ON THE WORK OF COLLE-SLAVETTI
C  [THEOR. CHIM. ACTA 37, 329(1975)]
C  IT IS DESCRIBED IN PHYSICAL REVIEW B 37, 785 (1988)
C  CHENGTEH LEE, WEITAO YANG, AND ROBERT G. PARR.
C  IT IS EXPLAINED MORE IN CHEM PHYS LETT 157, 200
C  (1989) MIEHLICH, SAVIN, STOLL, AND PREUSS.
C  POPLE, GILL AND JOHNSON FURTHER EXPLAIN IT IN
C  J CHEM PHYS 98, 5612 (1993).
C  THIS FUNCTIONAL DOESN'T DEPEND ON A "LOCAL" FUNCTIONAL AS ITS BASIS
C  THIS FUNCTIONAL HAS ITS BASIS IN THE SECOND ORDER DENSITY MATRIX
C
C  THIS FUNCTIONAL IS IMPLEMENTED FOR THE OPEN SHELL CASE IN A NON-SPIN
C  POLARIZED MANNER.
C
C
      IF(HISPIN)THEN
        CALL VADD(R,1,ROPEN,1,R,1,L3)
      ELSE
        CALL DSCAL(L3,2.0D+00,R,1)
      END IF
      CALL DFTJAC(R,UM,NUM,TM,IWRK,WRK,S)
      EIGCUT=1.0D-09
      DO 830 I=1,NUM
        IF (R(I,1).LT.EIGCUT) R(I,1) = EIGCUT
        IF (R(I,1).LT.EIGCUT) THEN
          R3   (I,1)=0.0D+00
          RSA  (I,1)=0.0D+00
          RSAVE(I,1)=0.0D+00
          RY2  (I,1)=0.0D+00
          S    (I,1)=0.0D+00
          R    (I,1)=0.0D+00
        ELSE
          RHO=R(I,1)
          DUM=RHO**(1.0D+00/3.0D+00)
          DUM2=1.0D+00/DUM
          R3(I,1)  = DUM2
C         THIS IS THE W(RHO) FUNCTION FROM LYP
          RHO21=((DUM2*EXP(PARAM1*DUM2))**11)/(1.0D+00+D*DUM2)
C         THIS IS THE D(RHO) FUNCTION FROM LYP
          RHO23=(PARAM2+PARAM3*DUM2)/(DUM+D)
C         THIS THE DERIVATIVE OF THE W(RHO) FUNCTION FROM LYP
          RHO22=((DUM2*EXP(PARAM4*DUM2))**15)*
     *          (PARAM7+PARAM5*DUM2+PARAM6*DUM2**2)/
     *          (DUM2*(1.0D+00+D*DUM2)**2)
C         THIS IS THE DERIVATIVE OF THE D(RHO) FUNCTION FROM LYP
          RHO24=(PARAM8*DUM2**5+PARAM9*DUM2**4+PARAM10*DUM2**6)/
     *          (1.0D+00+D*DUM2)**2
          RSA  (I,1)=-A*RHO/(RHO+D*DUM**2)-A*B*RHO21*CF*DUM**11
C         ADD IN THE RHO**(8.0D+00/3.0D+00) SINCE WE USE X NOT GRAD(RHO)
          RSAVE(I,1)=-A*B*RHO21*RHO*DUM**8*
     *               (-1.0D+00/24.0D+00-7.0D+00/72.0D+00*RHO23)
C         ADD IN RHO**(4.0D+00/3.0D+00) SINCE WE USE X AND NOT GRAD(RHO)
          RY2(I,1)=(A*B/18.0D+00)*RHO21*(3.0D+00+7.0D+00*RHO23)*DUM**10
C         ONCE AGAIN RHO^(8/3) BECAUSE WE USED X NOT GRADIENT(RHO)
          S(I,1)=(A*B/144.0D+00*(RHO22*RHO*RHO*(6.0D+00+14.0D+00*
     *             RHO23) + RHO21*RHO*(12.0D+00+28.0D+00*
     *             RHO23+14.0D+00*RHO*RHO24)))*DUM**8
C  NOW THE (D LYP/D RHO) TERM WITH NO GRADIENT DEPENDENCE AT ALL
          R(I,1)=-A/(1.0D+00+D/DUM)*
     *           (D/(3.0D+00*(DUM+D))+1.0D+00) -
     *           ABCF*(RHO22*(DUM**14)*
     *           (0.5D+00**(11.0D+00/3.0D+00))+RHO21*
     *           (DUM**11)*(7.0D+00/3.0D+00)*
     *           (0.5D+00)**(8.0D+00/3.0D+00))
        END IF
 830  CONTINUE
C
        CALL DFTTRN1(R3,UM,TM,NUM)
        CALL DFTTRN1(RSA,UM,TM,NUM)
        CALL DFTTRN1(RSAVE,UM,TM,NUM)
        CALL DFTTRN1(RY2,UM,TM,NUM)
        CALL DFTTRN1(S,UM,TM,NUM)
        CALL DFTTRN1(R,UM,TM,NUM)
C
C     RSA IS THE NON-GRADIENT TERMS
C     RSAVE IS THE GRADIENT DEPENDENT TERMS
C
C     GENERATE X
C
        DO 201 J=1,NUM
        DO 201 I=1,J
        RTX(I,J)=0.0D+00
        RTY(I,J)=0.0D+00
        RTZ(I,J)=0.0D+00
        DO 202 M=1,NUM
C
C      TAKE ADVANTAGE OF ANTI-SYMMETRIC DX AND SYMMETRIC R3 TO SPEED
C      UP MEMORY ACCESS
C       RTX(I,J)=RTX(I,J) -3.0D+00*(DX(I,M)*R3(M,J)-R3(I,M)*DX(M,J))
C
        RTX(I,J)=RTX(I,J) +(DX(M,I)*R3(M,J)+R3(M,I)*DX(M,J))
        RTY(I,J)=RTY(I,J) +(DY(M,I)*R3(M,J)+R3(M,I)*DY(M,J))
        RTZ(I,J)=RTZ(I,J) +(DZ(M,I)*R3(M,J)+R3(M,I)*DZ(M,J))
 202    CONTINUE
        RTX(I,J)=RTX(I,J)*3.0D+00
        RTY(I,J)=RTY(I,J)*3.0D+00
        RTZ(I,J)=RTZ(I,J)*3.0D+00
        RTX(J,I)=RTX(I,J)
        RTY(J,I)=RTY(I,J)
        RTZ(J,I)=RTZ(I,J)
 201    CONTINUE
C
C   RTZ=M[(D RHO/DZ)/(RHO**(4/3))]
C
      DO I=1,NUM
        DO J=1,I
          RX2(I,J)=0.0D+00
          DO M=1,NUM
            RX2(I,J)=RX2(I,J)+RTX(M,I)*RTX(M,J)
     *                       +RTY(M,I)*RTY(M,J)
     *                       +RTZ(M,I)*RTZ(M,J)
          END DO
          RX2(J,I)=RX2(I,J)
        END DO
      END DO
C
C   RX2=|GRADIENT RHO/RHO^4/3|^2
C
C   VERY IMPORTANT.  YOU MUST USE "X" RATHER THAN JUST THE GRADIENT SO
C   THAT THE GRADIENT TERMS ARE DIMENSIONLESS.  YOU LATER MULTIPLY BY
C   RHO^4/3 OR RHO^8/3 TO GET THE CORRECT EQUATION.
C
        DO 316 J=1,NUM
        DO 316 I=1,NUM
        RZ2(I,J)=0.0D+00
        DO 316 N=1,NUM
C        RZ2(I,J)=RZ2(I,J)+RTX(I,N)*DX(N,J)+RTY(I,N)*DY(N,J)+
C     *           RTZ(I,N)*DZ(N,J)
C        REVERSE RTX INDEX TO SPEED UP MEMORY ACCESS
         RZ2(I,J)=RZ2(I,J)+RTX(N,I)*DX(N,J)
     *                    +RTY(N,I)*DY(N,J)
     *                    +RTZ(N,I)*DZ(N,J)
 316    CONTINUE
C
C   RZ2(I,J)=(I (1/RHO^4/3) GRADIENT RHO GRADIENT J)
C
        DO 300 I=1,NUM
         DO 301 J=1,I
          TEMP=0.0D+00
          DO 302 M=1,NUM
C          TEMP=TEMP+(RSAVE(I,M)*RX2(M,J)+RX2(I,M)*RSAVE(M,J))/2.0D+00
           TEMP=TEMP+(RSAVE(M,I)*RX2(M,J)+RX2(M,I)*RSAVE(M,J))
 302      CONTINUE
          RSA(J,I)=RSA(J,I)+TEMP*0.5D+00
          RSA(I,J)=RSA(J,I)
 301     CONTINUE
 300    CONTINUE
C
C   NOW THE CLOSED SHELL FOCK MATRIX
C
C   START WITH THE DERIVATIVES WITH RESPECT TO THE GRADIENT
C
        DO 650 I=1,NUM
         DO 651 J=1,I
          RTX(J,I)=0.0D+00
          DO 652 M=1,NUM
C         NOTE THE RZ2 IS NOT SYMMETRIC BUT RY2 IS
          RTX(J,I)=RTX(J,I)+(RY2(M,I)*RZ2(M,J)+RY2(M,J)*RZ2(M,I))
 652      CONTINUE
C         MULTIPLY BY ONE-HALF IS DONE BELOW, SO IS SYMMETRIZING
C         RTX(J,I)=RTX(J,I)*0.5D+00
C         RTX(I,J)=RTX(J,I)
 651     CONTINUE
 650    CONTINUE
C
C   RTX IS NOW THE FINISHED FOCK FOR ALL TERMS THAT CONTAIN
C   DERIVATIVES WITH RESPECT TO THE GRADIENT
C   NOW WE DO THE DERIVATIVE WITH RESPECT TO DENSITY.
C   IT CONTAINS SECOND DERIVATIVE TERMS (EQ. A28 OF THE POPLE PAPER)
C
        DO 800 I=1,NUM
         DO 800 J=1,I
          DO 810 M=1,NUM
           RTX(J,I)=RTX(J,I)+(S(M,J)*RX2(M,I)+S(M,I)*RX2(M,J))
 810      CONTINUE
          R(J,I)=R(J,I)+RTX(J,I)*0.5D+00
          R(I,J)=R(J,I)
 800    CONTINUE
      IF(HISPIN)THEN
        CALL DCOPY(L3,RSA,1,RSB,1)
        CALL DCOPY(L3,R,1,ROPEN,1)
      END IF
      END
C*MODULE DFTFUN  *DECK DFTJAC
      SUBROUTINE DFTJAC(F,V,N,EIG,IWRK,WRK,S)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION F(N,N),V(N,N),EIG(N),IWRK(N),WRK(N,8),S(*)
C
C     F   IS A SQUARE AND SYMMETRIC THE MATRIX TO BE DIAGONALIZED
C     V   IS THE EIGENVECTORS TRANSPOSED
C     EIG IS THE EIGENVALUES
C     WRK,IWRK, AND S ARE ALL WORK AREAS
C
      CALL CPYSQT(F,S,N,1)
      CALL GLDIAG(N,N,N,S,WRK,EIG,V,IERR,IWRK)
      IF (IERR.NE.0) CALL ABRT
      DO 20 I=1,N
   20   F(I,1) = EIG(I)
C     TRANSPOSE EIGENVECTORS NOW
      CALL TRPOSQ(V,N)
      END
C*MODULE DFTFUN  *DECK DFTTRN1
      SUBROUTINE DFTTRN1(A,B,C,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(N,N),B(N,N),C(N)
C
C   THIS ROUTINE DOES A=(B DAGGER)*A*B ON DIAGONAL MATRIX
C   THIS ASSUMES THAT THE DIAGONAL VALUES ARE STORED IN (1,1) TO (N,1)
C   DFTJAC PRE-TRANSPOSED B, SO THAT WE DON'T DO B*A*(B DAGGER)
C   WHICH IS SLOWER
C
      DO I=1,N
        C(I)=A(I,1)
      END DO
      DO 200 I=1,N
        DO 200 J=1,I
          A(J,I)=0.0D+00
          DO 300 K=1,N
            A(J,I)=A(J,I)+B(K,I)*B(K,J)*C(K)
 300      CONTINUE
          A(I,J)=A(J,I)
 200  CONTINUE
      END
