C 14 NOV 05 - DGF - PAD COMMON BLOCK SCFOPT
C 19 SEP 05 - GDF - PAD SHLNRM COMMON FOR H AND I FUNCTIONS
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL
C 16 JAN 04 - DGF - NEW DFT GRADIENT THRESHOLDS
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 26 MAR 02 - DGF - ATMVEC: FIX FOR DUMMY ATOMS
C 24 JAN 02 - RMO - GRDPT,DFTTRF: FIX FOR AXP PROCESSOR
C  5 DEC 01 - DGF - FIX FOR CASE OF NO BETA ELECTRONS
C 25 JUN 01 - MWS - ALTER COMMON BLOCK SCFOPT AND WFNOPT
C 13 JUN 01 - DGF,MK ADD WEIGHT DERIVATIVES TO COMPLETE THE DFT GRADIENT
C 30 APR 01 - TT,SY,MK,DGF IMPLEMENT GRID-BASED DFT
C
C   POSSIBLE OPTIMIZATION: SOME ARRAYS LIKE ATMXVEC ARE ADDRESSED
C                          AS COLUMNS RATHER THAN ROWS.
C
C*MODULE DFTGRD  *DECK ATMVEC
      SUBROUTINE ATMVEC(ATMXVEC,ATMYVEC,ATMZVEC,RIJ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000,ZERO=0.0D+00)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      DIMENSION ATMXVEC(NAT,NAT),ATMYVEC(NAT,NAT),ATMZVEC(NAT,NAT),
     *          RIJ(NAT,NAT)
C
      DO 10 IATM=1,NAT
         XCNTR=C(1,IATM)
         YCNTR=C(2,IATM)
         ZCNTR=C(3,IATM)
         ATMXVEC(IATM,IATM)=ZERO
         ATMYVEC(IATM,IATM)=ZERO
         ATMZVEC(IATM,IATM)=ZERO
         RIJ(IATM,IATM)=ZERO
         DO 20 JATM=1,IATM-1
            ATMXVEC(IATM,JATM)= XCNTR-C(1,JATM)
            ATMYVEC(IATM,JATM)= YCNTR-C(2,JATM)
            ATMZVEC(IATM,JATM)= ZCNTR-C(3,JATM)
            ATMXVEC(JATM,IATM)=-ATMXVEC(IATM,JATM)
            ATMYVEC(JATM,IATM)=-ATMYVEC(IATM,JATM)
            ATMZVEC(JATM,IATM)=-ATMZVEC(IATM,JATM)
            RIJ(IATM,JATM)=SQRT(
     *            ATMXVEC(IATM,JATM)*ATMXVEC(IATM,JATM)
     *           +ATMYVEC(IATM,JATM)*ATMYVEC(IATM,JATM)
     *           +ATMZVEC(IATM,JATM)*ATMZVEC(IATM,JATM))
            RIJ(JATM,IATM)=RIJ(IATM,JATM)
 20      CONTINUE
 10   CONTINUE
      RETURN
      END
C*MODULE DFTGRD  *DECK DFTAO
      SUBROUTINE DFTAO(IJX,IJY,IJZ,VALMOA,ANGXVL,ANGYVL,ANGZVL,CL,EXPS,
     *                 RSQRD,NAT,L1)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, TWO=2.0D+00,RLN10=2.30258D+00)
C
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SHLNRM/ PNRM(84)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
C
      DIMENSION VALMOA(*),EXPS(*),ANGXVL(NAT,*),ANGYVL(NAT,*),
     >          ANGZVL(NAT,*),CL(*),RSQRD(NAT)
      DIMENSION IJX(35),IJY(35),IJZ(35)
C
C     THIS SUBROUTINE COMPUTES VALUES OF ALL AOS FOR A GRID POINT
C
C
      TOL = RLN10*ITOL/TWO
C
C     ITOL SETS THRESHOLD FOR INTEGRALS. HERE IT IS APPLIED TO
C     WAVEFUNCTION, THUS IT IS HALVED. IT IS FOUND THAT THIS IS A
C     RATHER CONSERVATIVE SETTING AND ONE GETS ABOUT 11 SIG. FIGURES
C     OF ACCURACY IN THE TOTAL ENERGY WITH THE DEFAULT ITOL=20.
C     SEE DFTGAO AND DFTG2AO.
C
      IFCT=0
      CALL VCLR(VALMOA,1,L1)
      DO 40 ISHELL = 1, NSHELL
         IATM = KATOM(ISHELL)
C***********************************************************************
C     ANGULAR INTERMEDIATES FOR THE DENSITY AND GRADIENT
C***********************************************************************
         K1 = KSTART(ISHELL)
         K2 = K1 + KNG(ISHELL) - 1
         MINI=KMIN(ISHELL)
         MAXI=KMAX(ISHELL)
         LOCI = KLOC(ISHELL)-MINI
         LOCI0= KLOC(ISHELL)
         DO 60 IMOMFCT = K1, K2
            IFCT=IFCT+1
            DUM=EX(IMOMFCT)*RSQRD(IATM)
            IF(DUM.GT.TOL) GOTO 60
            VEXP=EXP(-DUM)
            EXPS(IFCT)=VEXP
C           SAVE EXP FOR THE GRADIENT
            CC=CL(IMOMFCT)
            LOCI1=LOCI0
            IF(MINI.EQ.1) THEN
C
C           SPECIAL FAST CODE FOR S FUNCTIONS (INCLUDING L SHELLS)
C
               VALMOA(LOCI1)=VALMOA(LOCI1)+CC*VEXP
C              FOR L-SHELLS
               LOCI1=LOCI1+1
               CC=CP(IMOMFCT)
            ENDIF
            IF(MAXI.EQ.4) THEN
C
C           SPECIAL FAST CODE FOR P FUNCTIONS (INCLUDING L SHELLS)
C
               VALMOA(LOCI1  )=VALMOA(LOCI1  )+CC*VEXP*ANGXVL(IATM,3)
               VALMOA(LOCI1+1)=VALMOA(LOCI1+1)+CC*VEXP*ANGYVL(IATM,3)
               VALMOA(LOCI1+2)=VALMOA(LOCI1+2)+CC*VEXP*ANGZVL(IATM,3)
               LOCI1=LOCI1+3
            ENDIF
            IF(LOCI1.NE.LOCI0) GOTO 60
C
C           GENERAL CODE (WORKS FOR S,P,L FUNCTIONS AS WELL)
C           TO USE FOR L-FUNCTIONS, UNCOMMENT THE "C     L-SHELL" LINE
C
            DO 50 ITYP = MINI,MAXI
               IX=IJX(ITYP)-1
               IY=IJY(ITYP)-1
               IZ=IJZ(ITYP)-1
C
C              COMPUTE AO VALUE AT A GRID POINT
C
               DUM=CC*PNRM(ITYP)*VEXP
               VALMOA(LOCI+ITYP)=VALMOA(LOCI+ITYP)+ANGXVL(IATM,IX+2)*
     *            ANGYVL(IATM,IY+2)*ANGZVL(IATM,IZ+2)*DUM
C              UNCOMMENT FOR L-SHELL
C              IF(ITYP.EQ.1) CC=CP(IMOMFCT)
C              THIS IS THE DREADED L-SHELL CASE, FOR S-SHELLS THE LOOP
C              IS OVER, OTHERWISE USE P COEFFICIENT IN THE FOLLOWING
 50         CONTINUE
 60      CONTINUE
 40   CONTINUE
C
      RETURN
      END
C*MODULE DFTGRD  *DECK DFTGAO
      SUBROUTINE DFTGAO(IJX,IJY,IJZ,GVALAX,GVALAY,GVALAZ,ANGXVL,ANGYVL,
     *                  ANGZVL,CL,EXPS,RSQRD,NAT,L1)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, TWO=2.0D+00,RLN10=2.30258D+00)
C
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SHLNRM/ PNRM(84)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
C
      DIMENSION GVALAX(L1),GVALAY(L1),GVALAZ(L1),ANGXVL(NAT,*),
     >          ANGYVL(NAT,*),ANGZVL(NAT,*),CL(*),RSQRD(NAT),EXPS(*)
      DIMENSION IJX(35),IJY(35),IJZ(35)
C
C     THIS SUBROUTINE COMPUTES VALUES OF GRADIENT (DERIVATIVE BY THE
C     ELECTRON COORDINATE) OF ALL AOS FOR A GRID POINT
C
C
      TOL = RLN10*ITOL/TWO
C     THE SAME AS IN DFTAO.
C
      IFCT=0
      CALL VCLR(GVALAX,1,L1)
      CALL VCLR(GVALAY,1,L1)
      CALL VCLR(GVALAZ,1,L1)
      DO 40 ISHELL = 1, NSHELL
         IATM = KATOM(ISHELL)
C***********************************************************************
C     ANGULAR INTERMEDIATES FOR THE DENSITY AND GRADIENT
C***********************************************************************
         K1 = KSTART(ISHELL)
         K2 = K1 + KNG(ISHELL) - 1
         MINI=KMIN(ISHELL)
         MAXI=KMAX(ISHELL)
         LOCI = KLOC(ISHELL)-MINI
         LOCI0= KLOC(ISHELL)
         DO 60 IMOMFCT = K1, K2
            IFCT=IFCT+1
            TWOALP = EX(IMOMFCT)*TWO
            DUM=EX(IMOMFCT)*RSQRD(IATM)
            IF(DUM.GT.TOL) GOTO 60
C           VEXP=EXP(-DUM)
            VEXP=EXPS(IFCT)
            CC=CL(IMOMFCT)
            LOCI1=LOCI0
            IF(MINI.EQ.1) THEN
C
C           SPECIAL FAST CODE FOR S FUNCTIONS (INCLUDING L SHELLS)
C
              DUM=CC*VEXP*TWOALP
              GVALAX(LOCI1)=GVALAX(LOCI1)-DUM*ANGXVL(IATM,3)
              GVALAY(LOCI1)=GVALAY(LOCI1)-DUM*ANGYVL(IATM,3)
              GVALAZ(LOCI1)=GVALAZ(LOCI1)-DUM*ANGZVL(IATM,3)
C             FOR L-SHELLS
              LOCI1=LOCI1+1
              CC=CP(IMOMFCT)
            ENDIF
            IF(MAXI.EQ.4) THEN
C
C           SPECIAL FAST CODE FOR P FUNCTIONS (INCLUDING L SHELLS)
C
               DUM=CC*TWOALP*VEXP
               DUM1=CC*VEXP
               XY=DUM*ANGXVL(IATM,3)*ANGYVL(IATM,3)
               XZ=DUM*ANGXVL(IATM,3)*ANGZVL(IATM,3)
               YZ=DUM*ANGYVL(IATM,3)*ANGZVL(IATM,3)
               GVALAX(LOCI1)=GVALAX(LOCI1)+DUM1-DUM*ANGXVL(IATM,4)
               GVALAY(LOCI1)=GVALAY(LOCI1)-XY
               GVALAZ(LOCI1)=GVALAZ(LOCI1)-XZ
               LOCI1=LOCI1+1
               GVALAX(LOCI1)=GVALAX(LOCI1)-XY
               GVALAY(LOCI1)=GVALAY(LOCI1)+DUM1-DUM*ANGYVL(IATM,4)
               GVALAZ(LOCI1)=GVALAZ(LOCI1)-YZ
               LOCI1=LOCI1+1
               GVALAX(LOCI1)=GVALAX(LOCI1)-XZ
               GVALAY(LOCI1)=GVALAY(LOCI1)-YZ
               GVALAZ(LOCI1)=GVALAZ(LOCI1)+DUM1-DUM*ANGZVL(IATM,4)
               LOCI1=LOCI1+1
            ENDIF
            IF(LOCI1.NE.LOCI0) GOTO 60
C
C           GENERAL CODE (WORKS FOR S,P,L FUNCTIONS AS WELL)
C           TO USE FOR L-FUNCTIONS, UNCOMMENT THE "C     L-SHELL" LINE
C
            DO 50 ITYP = MINI,MAXI
               IX=IJX(ITYP)-1
               IY=IJY(ITYP)-1
               IZ=IJZ(ITYP)-1
               DUM=CC*PNRM(ITYP)*VEXP
C
C              COMPUTE GRADIENT AO VALUE AT A GRID POINT
C              GRADIENT IS BY THE ELECTRON (NOT NUCLEAR) COORDINATES
C
               ANGYZVL=ANGYVL(IATM,IY+2)*ANGZVL(IATM,IZ+2)
               ANGZXVL=ANGZVL(IATM,IZ+2)*ANGXVL(IATM,IX+2)
               ANGXYVL=ANGXVL(IATM,IX+2)*ANGYVL(IATM,IY+2)
C***********************************************************************
C     GRADIENT MINUS ONE COMPONENT
C***********************************************************************
               ANGXM= IX*ANGXVL(IATM,IX+1)*ANGYZVL
               ANGYM= IY*ANGYVL(IATM,IY+1)*ANGZXVL
               ANGZM= IZ*ANGZVL(IATM,IZ+1)*ANGXYVL
C***********************************************************************
C     GRADIENT PLUS ONE COMPONENT
C***********************************************************************
               ANGXP= ANGXVL(IATM,IX+3)*ANGYZVL
               ANGYP= ANGYVL(IATM,IY+3)*ANGZXVL
               ANGZP= ANGZVL(IATM,IZ+3)*ANGXYVL
C
               LOCI0=LOCI+ITYP
               GVALAX(LOCI0) = GVALAX(LOCI0)+DUM*(ANGXM-TWOALP*ANGXP)
               GVALAY(LOCI0) = GVALAY(LOCI0)+DUM*(ANGYM-TWOALP*ANGYP)
               GVALAZ(LOCI0) = GVALAZ(LOCI0)+DUM*(ANGZM-TWOALP*ANGZP)
C              UNCOMMENT FOR L-SHELLS
C              IF(ITYP.EQ.1) CC=CP(IMOMFCT)
C              THIS IS THE DREADED L-SHELL CASE, FOR S-SHELLS THE LOOP
C              IS OVER, OTHERWISE USE P COEFFICIENT IN THE FOLLOWING.
 50         CONTINUE
 60      CONTINUE
 40   CONTINUE
C
      RETURN
      END
C*MODULE DFTGRD  *DECK DFTG2AO
      SUBROUTINE DFTG2AO(IJX,IJY,IJZ,G2AOXX,G2AOYY,G2AOZZ,G2AOXY,G2AOYZ,
     *                 G2AOXZ,ANGXVL,ANGYVL,ANGZVL,CL,EXPS,RSQRD,NAT,L1)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, ONE=1.0D+00,
     >             TWO=2.0D+00,THREE=3.0D+00,RLN10=2.30258D+00)
C
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SHLNRM/ PNRM(84)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
C
      DIMENSION CL(*),RSQRD(NAT),G2AOXX(*),G2AOYY(*),G2AOZZ(*),
     >          G2AOXY(*),G2AOYZ(*),G2AOXZ(*),EXPS(*),
     >          ANGXVL(NAT,*),ANGYVL(NAT,*),ANGZVL(NAT,*)
      DIMENSION IJX(35),IJY(35),IJZ(35)
C
C     COMPUTE SECOND DERIVATIVES OF THE BASIS FUNCTIONS
C
      TOL = RLN10*ITOL/TWO
C     THE SAME AS IN DFTAO.
C
      IFCT=0
      CALL VCLR(G2AOXX,1,L1)
      CALL VCLR(G2AOYY,1,L1)
      CALL VCLR(G2AOZZ,1,L1)
      CALL VCLR(G2AOXY,1,L1)
      CALL VCLR(G2AOYZ,1,L1)
      CALL VCLR(G2AOXZ,1,L1)
      DO 40 ISHELL = 1, NSHELL
         IATM = KATOM(ISHELL)
C
C     ANGULAR INTERMEDIATES FOR THE DENSITY AND GRADIENT
C
         K1 = KSTART(ISHELL)
         K2 = K1 + KNG(ISHELL) - 1
         MINI=KMIN(ISHELL)
         MAXI=KMAX(ISHELL)
         LOCI = KLOC(ISHELL)-MINI
         LOCI0= KLOC(ISHELL)
         DO 60 IMOMFCT = K1, K2
           IFCT=IFCT+1
           TWOALP = EX(IMOMFCT)*TWO
           DUM=EX(IMOMFCT)*RSQRD(IATM)
           IF(DUM.GT.TOL) GOTO 60
           VEXP=EXPS(IFCT)
           CC=CL(IMOMFCT)
           LOCI1=LOCI0
           IF(MINI.EQ.1) THEN
C
C     SPECIAL FAST CODE FOR S FUNCTIONS (INCLUDING L SHELLS)
C
             DUM=CC*VEXP*TWOALP
C
             XY=DUM*ANGXVL(IATM,3)*ANGYVL(IATM,3)
             XZ=DUM*ANGXVL(IATM,3)*ANGZVL(IATM,3)
             YZ=DUM*ANGYVL(IATM,3)*ANGZVL(IATM,3)
C
             G2AOXX(LOCI1)=G2AOXX(LOCI1)-DUM*(ONE-TWOALP*ANGXVL(IATM,4))
             G2AOYY(LOCI1)=G2AOYY(LOCI1)-DUM*(ONE-TWOALP*ANGYVL(IATM,4))
             G2AOZZ(LOCI1)=G2AOZZ(LOCI1)-DUM*(ONE-TWOALP*ANGZVL(IATM,4))
             G2AOXY(LOCI1)=G2AOXY(LOCI1)+XY*TWOALP
             G2AOYZ(LOCI1)=G2AOYZ(LOCI1)+YZ*TWOALP
             G2AOXZ(LOCI1)=G2AOXZ(LOCI1)+XZ*TWOALP
C
C     FOR L-SHELLS
C
             LOCI1=LOCI1+1
             CC=CP(IMOMFCT)
           ENDIF
           IF(MAXI.EQ.4) THEN
C
C     SPECIAL FAST CODE FOR P FUNCTIONS (INCLUDING L SHELLS)
C
             DUM=CC*TWOALP*VEXP
             DUMX=ONE-TWOALP*ANGXVL(IATM,4)
             DUMY=ONE-TWOALP*ANGYVL(IATM,4)
             DUMZ=ONE-TWOALP*ANGZVL(IATM,4)
             DUMXYZ=TWOALP*DUM*ANGXVL(IATM,3)*ANGYVL(IATM,3)*
     *                         ANGZVL(IATM,3)
             TWOALP2=TWOALP*TWOALP
C
C     PX
C     ----
             G2AOXX(LOCI1)=G2AOXX(LOCI1)
     >                 -DUM*(THREE*ANGXVL(IATM,3)-TWOALP*ANGXVL(IATM,5))
             G2AOYY(LOCI1)=G2AOYY(LOCI1)-DUM*DUMY*ANGXVL(IATM,3)
             G2AOZZ(LOCI1)=G2AOZZ(LOCI1)-DUM*DUMZ*ANGXVL(IATM,3)
C
             G2AOXY(LOCI1)=G2AOXY(LOCI1)-DUM*DUMX*ANGYVL(IATM,3)
             G2AOYZ(LOCI1)=G2AOYZ(LOCI1)+DUMXYZ
             G2AOXZ(LOCI1)=G2AOXZ(LOCI1)-DUM*DUMX*ANGZVL(IATM,3)
             LOCI1=LOCI1+1
C
C     PY
C     ----
             G2AOXX(LOCI1)=G2AOXX(LOCI1)-DUM*DUMX*ANGYVL(IATM,3)
             G2AOYY(LOCI1)=G2AOYY(LOCI1)
     >                 -DUM*(THREE*ANGYVL(IATM,3)-TWOALP*ANGYVL(IATM,5))
             G2AOZZ(LOCI1)=G2AOZZ(LOCI1)-DUM*DUMZ*ANGYVL(IATM,3)
C
             G2AOXY(LOCI1)=G2AOXY(LOCI1)-DUM*DUMY*ANGXVL(IATM,3)
             G2AOYZ(LOCI1)=G2AOYZ(LOCI1)-DUM*DUMY*ANGZVL(IATM,3)
             G2AOXZ(LOCI1)=G2AOXZ(LOCI1)+DUMXYZ
             LOCI1=LOCI1+1
C
C     PZ
C     ----
             G2AOXX(LOCI1)=G2AOXX(LOCI1)-DUM*DUMX*ANGZVL(IATM,3)
             G2AOYY(LOCI1)=G2AOYY(LOCI1)-DUM*DUMY*ANGZVL(IATM,3)
             G2AOZZ(LOCI1)=G2AOZZ(LOCI1)
     >                 -DUM*(THREE*ANGZVL(IATM,3)-TWOALP*ANGZVL(IATM,5))
C
             G2AOXY(LOCI1)=G2AOXY(LOCI1)+DUMXYZ
             G2AOYZ(LOCI1)=G2AOYZ(LOCI1)-DUM*DUMZ*ANGYVL(IATM,3)
             G2AOXZ(LOCI1)=G2AOXZ(LOCI1)-DUM*DUMZ*ANGXVL(IATM,3)
             LOCI1=LOCI1+1
C
           ENDIF
           IF(LOCI1.NE.LOCI0) GOTO 60
C
C     GENERAL CODE (WORKS FOR S,P,L FUNCTIONS AS WELL)
C     TO USE FOR L-FUNCTIONS, UNCOMMENT THE "C     L-SHELL" LINE
C
           DO 50 ITYP = MINI,MAXI
             IX=IJX(ITYP)-1
             IY=IJY(ITYP)-1
             IZ=IJZ(ITYP)-1
             DUM=CC*PNRM(ITYP)*VEXP
C
C            COMPUTE 2ND DERIVATIVE OF AO VALUE AT A GRID POINT
C            GRADIENT IS BY THE ELECTRON (NOT NUCLEAR) COORDINATES
C
             ANGYZVL=ANGYVL(IATM,IY+2)*ANGZVL(IATM,IZ+2)
             ANGZXVL=ANGZVL(IATM,IZ+2)*ANGXVL(IATM,IX+2)
             ANGXYVL=ANGXVL(IATM,IX+2)*ANGYVL(IATM,IY+2)
C
C     GRADIENT MINUS TWO COMPONENT
C     ---------------------------
             ANGXMM= IX*(IX-1)*ANGXVL(IATM,IX)*ANGYZVL
             ANGYMM= IY*(IY-1)*ANGYVL(IATM,IY)*ANGZXVL
             ANGZMM= IZ*(IZ-1)*ANGZVL(IATM,IZ)*ANGXYVL
             ANGXYM= IX*IY*ANGXVL(IATM,IX+1)
     >                    *ANGYVL(IATM,IY+1)*ANGZVL(IATM,IZ+2)
             ANGXZM= IX*IZ*ANGXVL(IATM,IX+1)
     >                    *ANGYVL(IATM,IY+2)*ANGZVL(IATM,IZ+1)
             ANGYZM= IY*IZ*ANGXVL(IATM,IX+2)
     >                    *ANGYVL(IATM,IY+1)*ANGZVL(IATM,IZ+1)
C
C     GRADIENT PLUS TWO COMPONENT
C     --------------------------
             ANGXPP= ANGXVL(IATM,IX+4)*ANGYZVL
             ANGYPP= ANGYVL(IATM,IY+4)*ANGZXVL
             ANGZPP= ANGZVL(IATM,IZ+4)*ANGXYVL
           ANGXYP= ANGXVL(IATM,IX+3)*ANGYVL(IATM,IY+3)*ANGZVL(IATM,IZ+2)
           ANGXZP= ANGXVL(IATM,IX+3)*ANGYVL(IATM,IY+2)*ANGZVL(IATM,IZ+3)
           ANGYZP= ANGXVL(IATM,IX+2)*ANGYVL(IATM,IY+3)*ANGZVL(IATM,IZ+3)
C
C     OTHERS
C     ------
          ANGXY=IY*ANGXVL(IATM,IX+3)*ANGYVL(IATM,IY+1)*ANGZVL(IATM,IZ+2)
          ANGYX=IX*ANGXVL(IATM,IX+1)*ANGYVL(IATM,IY+3)*ANGZVL(IATM,IZ+2)
          ANGXZ=IZ*ANGXVL(IATM,IX+3)*ANGYVL(IATM,IY+2)*ANGZVL(IATM,IZ+1)
          ANGZX=IX*ANGXVL(IATM,IX+1)*ANGYVL(IATM,IY+2)*ANGZVL(IATM,IZ+3)
          ANGYZ=IZ*ANGXVL(IATM,IX+2)*ANGYVL(IATM,IY+3)*ANGZVL(IATM,IZ+1)
          ANGZY=IY*ANGXVL(IATM,IX+2)*ANGYVL(IATM,IY+1)*ANGZVL(IATM,IZ+3)
C
C
          ANGDEN=ANGXVL(IATM,IX+2)*ANGYVL(IATM,IY+2)*ANGZVL(IATM,IZ+2)
C
          TWOALP2= TWOALP * TWOALP
C
          LOCI0=LOCI+ITYP
          G2AOXX(LOCI0)=G2AOXX(LOCI0)+DUM*(ANGXMM-TWOALP*(IX*2+1)*ANGDEN
     *                                   + TWOALP2*ANGXPP)
          G2AOYY(LOCI0)=G2AOYY(LOCI0)+DUM*(ANGYMM-TWOALP*(IY*2+1)*ANGDEN
     >                                   + TWOALP2*ANGYPP)
          G2AOZZ(LOCI0)=G2AOZZ(LOCI0)+DUM*(ANGZMM-TWOALP*(IZ*2+1)*ANGDEN
     >                                   + TWOALP2*ANGZPP)
          G2AOXY(LOCI0)=G2AOXY(LOCI0)+DUM*(ANGXYM-TWOALP*(ANGXY+ANGYX)
     >                                   + TWOALP2*ANGXYP)
          G2AOYZ(LOCI0)=G2AOYZ(LOCI0)+DUM*(ANGYZM-TWOALP*(ANGYZ+ANGZY)
     >                                   + TWOALP2*ANGYZP)
          G2AOXZ(LOCI0)=G2AOXZ(LOCI0)+DUM*(ANGXZM-TWOALP*(ANGXZ+ANGZX)
     >                                   + TWOALP2*ANGXZP)
C
C          L-SHELL
C          IF(ITYP.EQ.1) CC=CP(IMOMFCT)
C          THIS IS THE DREADED L-SHELL CASE, FOR S-SHELLS THE LOOP IS
C          OVER AND OTHERWISE USE THE P COEFFICIENT IN THE FOLLOWING.
 50        CONTINUE
 60      CONTINUE
 40   CONTINUE
      RETURN
      END
C*MODULE DFTGRD  *DECK DFTFOCK
      SUBROUTINE DFTFOCK(NEEDGR,UROHF,FTOTWT,DUMA,DUMB,DUMAX,DUMAY,
     *                   DUMAZ,DUMBX,DUMBY,DUMBZ,VALGA,VALGB,AOX,GAOX,
     *                   GAOY,GAOZ,FA,FB,CUTOFF,L1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF,NEEDGR
      PARAMETER (TWO=2.0D+00,ZERO=0.0D+00)
      DIMENSION AOX(*),GAOX(*),GAOY(*),GAOZ(*),FA(*),FB(*),VALGA(*),
     *          VALGB(*)
C
C     ADD DFT EXCHANGE/CORRELATION TO THE FOCK MATRIX (ALPHA AND BETA).
C
      IF(.NOT.NEEDGR) THEN
        IND=1
        DO I=1,L1
          FFA=FTOTWT*DUMA*AOX(I)
          IF(UROHF) THEN
            FFB=FTOTWT*DUMB*AOX(I)
            IF(ABS(FFA)+ABS(FFB).LT.CUTOFF*TWO) THEN
              IND=IND+I
              GOTO 100
            ENDIF
            DO J=1,I
              FA(IND)=FA(IND)+FFA*AOX(J)
              FB(IND)=FB(IND)+FFB*AOX(J)
              IND=IND+1
            ENDDO
  100       CONTINUE
C           CALL DAXPY(I,FFA,AOX,1,FA(IND),1)
C           CALL DAXPY(I,FFB,AOX,1,FB(IND),1)
C           IND=IND+I
          ELSE
            IF(ABS(FFA).LT.CUTOFF) THEN
              IND=IND+I
              GOTO 200
            ENDIF
            DO J=1,I
              FA(IND)=FA(IND)+FFA*AOX(J)
              IND=IND+1
            ENDDO
  200       CONTINUE
          ENDIF
        ENDDO
      ELSE
        IND=1
        AOMAX=ZERO
        AMAXVAL=ZERO
        BMAXVAL=ZERO
        DO I=1,L1
          VALGA(I)=(DUMA/TWO*AOX(I)+DUMAX*GAOX(I)+DUMAY*GAOY(I)+
     *              DUMAZ*GAOZ(I))*FTOTWT
          C1=AOX(I)
          C2A=VALGA(I)
          AOMAX=MAX(AOMAX,ABS(C1))
          AMAXVAL=MAX(AMAXVAL,ABS(C2A))
          IF(UROHF) THEN
           VALGB(I)=(DUMB/TWO*AOX(I)+DUMBX*GAOX(I)+DUMBY*GAOY(I)
     *              +DUMBZ*GAOZ(I))*FTOTWT
            C2B=VALGB(I)
            BMAXVAL=MAX(BMAXVAL,ABS(C2B))
            IF(ABS(C1*(AMAXVAL+BMAXVAL))+(ABS(C2A)+ABS(C2B))*AOMAX.LT.
     *         CUTOFF*TWO) THEN
              IND=IND+I
              GOTO 300
            ENDIF
            DO J=1,I
              FA(IND)=FA(IND)+C1*VALGA(J)+C2A*AOX(J)
              FB(IND)=FB(IND)+C1*VALGB(J)+C2B*AOX(J)
              IND=IND+1
            ENDDO
  300       CONTINUE
          ELSE
            IF(ABS(C1*AMAXVAL)+ABS(C2A*AOMAX).LT.CUTOFF) THEN
              IND=IND+I
              GOTO 400
            ENDIF
            DO J=1,I
              FA(IND)=FA(IND)+C1*VALGA(J)+C2A*AOX(J)
              IND=IND+1
            ENDDO
  400       CONTINUE
          ENDIF
        ENDDO
      ENDIF
      RETURN
      END
C*MODULE DFTGRD  *DECK DMATD
      SUBROUTINE DMATD(TOTWT,IIFACT,NAPTS,COEFFA,COEFFB,IANGN,IFACTR,
     *                 INC0,VMOA,DVMOXA,DVMOYA,DVMOZA,VALGA,VMOB,DVMOXB,
     *                 DVMOYB,DVMOZB,VALGB,FA,FB,NCNTR,EEXC,TOTELE,NANG,
     *                 PTRAD,XDAT,YDAT,ZDAT,ATMXVEC,ATMYVEC,ATMZVEC,
     *                 ANGXVL,ANGYVL,ANGZVL,BSLRDII,RSQRD,PCOEFF,EXPS,
     *                 AOX,GAOX,GAOY,GAOZ,IJX,IJY,IJZ,L1,BOHR,NEEDGR,
     *                 IZCORE,UROHF)
C***********************************************************************
C
C     THIS ROUTINE DOES A NUMERICAL INTEGRATION TO YIELD THE <I|V|A>
C     MATRIX WHERE V=D_E(XC)/D_RHO.  IN C1 SYMMETRY FOR NOW
C     THE RADIAL QUADRATURE FORMULA IS TAKEN FROM
C     P.M.W.GILL, B.G.JOHNSON, J.A.POPLE AND M.J.FRISCH,
C     CHEM. PHYS. LETT. 197, 499 (1992).
C     THE ANGULAR QUADRATURE FORMULA IS TAKEN FROM V.I.LEBEDEV,
C     ZH. VYCHISL. MAT. FIZ. 15, 48 (1975) AND 16, 293 (1976),
C     (ENGLISH TRANSLATION IN U.S.S.R. COMPUT. MATH AND MATH PHYS).
C
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL UROHF,OUT,NEEDGR,GOPARR,DSKWRK,MASWRK,NXT
C
      PARAMETER (MXATM=2000,ZERO=0.0D+00,ONE=1.0D+00,TWO=2.0D+00)
C
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
      DIMENSION TOTWT(NAT,*),COEFFA(L1,*),COEFFB(L1,*),VMOA(L1),
     *          DVMOXA(L1),DVMOYA(L1),DVMOZA(L1),VMOB(L1),DVMOXB(L1),
     *          DVMOYB(L1),DVMOZB(L1),VALGA(L1),VALGB(L1),IZCORE(*)
      DIMENSION NAPTS(NAT),IANGN(NAT,2),IIFACT(NAT),IFACTR(NAT)
      DIMENSION FA(*),FB(*),AOX(L1),GAOX(L1),GAOY(L1),GAOZ(L1),
     *          PTRAD(*),XDAT(NAT,NTHE*NPHI),YDAT(NAT,NTHE*NPHI),
     *          ZDAT(NAT,NTHE*NPHI),RSQRD(*),PCOEFF(*),BSLRDII(*),
     *          ATMXVEC(NAT,NAT),ATMYVEC(NAT,NAT),ATMZVEC(NAT,NAT),
     *          ANGXVL(NAT,*),ANGYVL(NAT,*),ANGZVL(NAT,*),EXPS(*)
      DIMENSION IJX(35),IJY(35),IJZ(35)
      DATA DEBUG/8HDEBUG   /, DFTGRD/8HDERDFT  /
C
C     EXCHANGE AND CORRELATION ENERGY CONTRIBUTION DUE TO INTEGRATION
C     OVER GRID POINTS SURROUNDING ATOM NCNTR IS COMPUTED.
C
      OUT = EXETYP.EQ.DFTGRD  .OR.  EXETYP.EQ.DEBUG
      NUC=INT(ZAN(NCNTR)+IZCORE(NCNTR))
      IF(NFG.NE.0) NUC=IAN(NCNTR)
      RAD=BSLRDII(NUC)/BOHR
      EXEC = ZERO
      NOA = NA
      NOB = NB
      NPT=NRAD*NTHE*NPHI
C
C     SET CUT-OFFS FOR THE DENSITY RCUTOFF AND WEIGHT WCUTOFF
C     RCUTOFF IS SET TO DEPEND UPON SCF DENSITY CONV. AND THE GRID SIZE
C     WCUTOFF IS A CELL VOLUME AND WE SET IT TO A FIXED VALUE.
C     MOST CELLS HAVE LARGE VOLUME (ABOUT 97% HAVE VOLUME .GT. 1E-14)
C
      DFTTHRS=DFTTHR
      IF(DFTTHR.EQ.ZERO) DFTTHR=1.0D-04/(NPT*NAT)
      WCUTOFF=1.0D-08/(NPT*NAT)
      RCUTOFF=CONVHF/(NPT*NAT)
      CCUTOFF=1.0D-03/(NPT*NAT)
      IF(DFTTHR.LT.1.1D-15) THEN
         WCUTOFF=1.0D-15
         RCUTOFF=1.0D-15
         CCUTOFF=1.0D-15
      ENDIF
C     WRITE(6,*) 'WWW',DFTTHR,WCUTOFF,RCUTOFF,CCUTOFF
C
C     ANGXVL=X**I, ANGYVL=Y**J, ANGZVL=Z**K, FOR NEEDED VALUES OF I,J,K
C     WHERE (X,Y,Z) IS THE CENTRE OF A DFT GRID POINT RI
C     MINUS AN ATOMIC CENTRE RA: XYZ= (RI - RA)
C
      CALL VCLR(ANGXVL(1,1),1,NAT)
      CALL VCLR(ANGYVL(1,1),1,NAT)
      CALL VCLR(ANGZVL(1,1),1,NAT)
      CALL DACOPY(NAT,ONE,ANGXVL(1,2),1)
      CALL DACOPY(NAT,ONE,ANGYVL(1,2),1)
      CALL DACOPY(NAT,ONE,ANGZVL(1,2),1)
C
C     ---- SYMMETRY MULTIPLICATION FACTOR ----
C
      FACT = IIFACT(NCNTR)*IFACTR(INC0)
C
C     DYNAMIC LOAD BALANCING: DIVIDE ALL POINTS INTO NGRAN*NPROC CHUNKS.
C     NGRAN CANNOT BE SET IN THE INPUT FILE. THEN EACH NODE GETS A CHUNK
C     AND TOILS IT OVER UNTIL READY TO ASK FOR MORE.
C
      LOOP=0
      NGRAN=20
      NLOOP=NRAD*(IANGN(NCNTR,2)-IANGN(NCNTR,1)+1)
      MCHUNK=(NLOOP-1)/(NGRAN*NPROC)+1
      NEXT  = -1
      NXT = IBTYP.EQ.1
C***********************************************************************
C     LOOP OVER RADIAL GRIDS
C     LOOP OVER ANGULAR GRIDS
C***********************************************************************
      DO 20 IRADPT = 1, NRAD
C
        R1= RAD*PTRAD(IRADPT)
        DO 10 IANGPT = IANGN(NCNTR,1), IANGN(NCNTR,2)
          IPT=(IRADPT-1)*NAPTS(NCNTR)+IANGPT
C
          IF(GOPARR) THEN
            IF(NXT) THEN
              LOOP=LOOP+1
              ICHUNK=(LOOP-1)/MCHUNK
              IF(ICHUNK.GT.NEXT) THEN
                CALL DDI_DLBNEXT(NEXT)
              ENDIF
              IF(NEXT.NE.ICHUNK) GOTO 10
            ELSE
              IF(MOD(IPT,NPROC).NE.ME) GOTO 10
            ENDIF
          ENDIF
          IPTME=(IPT-1)/NPROC+1
          IF(NXT) IPTME=IPT
C
          FTOTWT = TOTWT(NCNTR,IPTME)*FACT
          IF(FTOTWT.LT.WCUTOFF) THEN
             GOTO 10
          ENDIF
C
          XD=R1*XDAT(NCNTR,IANGPT)
          YD=R1*YDAT(NCNTR,IANGPT)
          ZD=R1*ZDAT(NCNTR,IANGPT)
          DO 610 IATM=1,NAT
            XCDNT=ATMXVEC(NCNTR,IATM)+XD
            YCDNT=ATMYVEC(NCNTR,IATM)+YD
            ZCDNT=ATMZVEC(NCNTR,IATM)+ZD
            RSQRD(IATM)=XCDNT**2+YCDNT**2+ZCDNT**2
            ANGXVL(IATM,3)=XCDNT
            ANGYVL(IATM,3)=YCDNT
            ANGZVL(IATM,3)=ZCDNT
            DO 35 IANG=3,NANG
               ANGXVL(IATM,IANG+1)=ANGXVL(IATM,IANG)*XCDNT
               ANGYVL(IATM,IANG+1)=ANGYVL(IATM,IANG)*YCDNT
               ANGZVL(IATM,IANG+1)=ANGZVL(IATM,IANG)*ZCDNT
 35         CONTINUE
  610     CONTINUE
          CALL DFTAO(IJX,IJY,IJZ,AOX,ANGXVL,ANGYVL,ANGZVL,PCOEFF,EXPS,
     *               RSQRD,NAT,L1)
C***********************************************************************
C     FORM DENSITY AT THIS POINT
C***********************************************************************
          CALL DFTTRF(UROHF,L1,NOA,NOB,COEFFA,COEFFB,AOX,VMOA,VMOB,
     *                CCUTOFF)
          ROA=DDOT(NOA,VMOA,1,VMOA,1)
          IF(UROHF) THEN
            ROB=DDOT(NOB,VMOB,1,VMOB,1)
          ELSE
            ROB=ROA
          ENDIF
          IF(ROA+ROB.LT.RCUTOFF) THEN
             GOTO 10
          ENDIF
          IF(NEEDGR) THEN
C***********************************************************************
C     FORM DENSITY GRADIENT AT THIS POINT
C***********************************************************************
            CALL DFTGAO(IJX,IJY,IJZ,GAOX,GAOY,GAOZ,ANGXVL,ANGYVL,ANGZVL,
     *                  PCOEFF,EXPS,RSQRD,NAT,L1)
            CALL DFTTRFG(UROHF,L1,NOA,NOB,COEFFA,COEFFB,
     *                   AOX,GAOX,GAOY,GAOZ,
     *                   DVMOXA,DVMOYA,DVMOZA,DVMOXB,DVMOYB,DVMOZB,
     *                   CCUTOFF)
            GRADXA=TWO*DDOT(NOA,VMOA,1,DVMOXA,1)
            GRADYA=TWO*DDOT(NOA,VMOA,1,DVMOYA,1)
            GRADZA=TWO*DDOT(NOA,VMOA,1,DVMOZA,1)
            IF(UROHF) THEN
              GRADXB=TWO*DDOT(NOB,VMOB,1,DVMOXB,1)
              GRADYB=TWO*DDOT(NOB,VMOB,1,DVMOYB,1)
              GRADZB=TWO*DDOT(NOB,VMOB,1,DVMOZB,1)
            ELSE
              GRADXB=GRADXA
              GRADYB=GRADYA
              GRADZB=GRADZA
            ENDIF
C***********************************************************************
C      FORM GRADIENT INVARIENT (GRAD DOT GRAD)
C***********************************************************************
            GRDAA  = GRADXA*GRADXA+GRADYA*GRADYA+GRADZA*GRADZA
            GRDBB  = GRADXB*GRADXB+GRADYB*GRADYB+GRADZB*GRADZB
            GRDAB  = GRADXA*GRADXB+GRADYA*GRADYB+GRADZA*GRADZB
          ENDIF
C***********************************************************************
C     THE EXCHANGE-CORRELATION FUNCTIONAL CALCULATION
C***********************************************************************
          VXCA1=ZERO
          VXCB1=ZERO
          DUMAX=ZERO
          DUMAY=ZERO
          DUMAZ=ZERO
          DUMBX=ZERO
          DUMBY=ZERO
          DUMBZ=ZERO
          XALPHA=ZERO
          XGRD=ZERO
          ECF=ZERO
          CALL CALCEXC(ROA,ROB,FTOTWT,GRDAA,GRDBB,GRDAB,
     >                 GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                 XALPHA,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >                 VXCB1,DUMBX,DUMBY,DUMBZ,ECF)
          DUMA = VXCA1
          DUMB = VXCB1
          EXEC1= XALPHA + XGRD + ECF
          EXEC = EXEC + EXEC1
          IF(NOB.EQ.0) THEN
             DUMB=ZERO
             GRADXB=ZERO
             GRADYB=ZERO
             GRADZB=ZERO
          ENDIF
C***********************************************************************
C      CONSTRUCT FOCK MATRIX
C***********************************************************************
          CALL DFTFOCK(NEEDGR,UROHF,FTOTWT,DUMA,DUMB,DUMAX,DUMAY,DUMAZ,
     *                 DUMBX,DUMBY,DUMBZ,VALGA,VALGB,AOX,GAOX,GAOY,GAOZ,
     *                 FA,FB,DFTTHR,L1)
C
C     ----- THE TOTAL ELECTRON DENSITY -----
C
          TOTELE=TOTELE+FTOTWT*(ROA+ROB)
C
   10   CONTINUE
C
C     ----- NEXT RADIAL POINT -----
C
   20 CONTINUE
C
C     ----- NEXT ATOM -----
C
      IF(GOPARR.AND.NXT) CALL DDI_DLBRESET
      EEXC = EEXC+EXEC
      DFTTHR=DFTTHRS
      IF (OUT) WRITE(IW,9999)NCNTR,EEXC,TOTELE
 9999 FORMAT(/5X,'ATM',I8,' EXC=',F20.10,5X,'TOTELE=',F20.10)
      RETURN
      END
C*MODULE DFTGRD  *DECK GLGQUD
      SUBROUTINE GLGQUD(X1,X2,X,W,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (EPS =3.0D-14)
      PARAMETER (ZERO=0.0D+00)
      PARAMETER (ONE =1.0D+00)
      PARAMETER (TWO =2.0D+00)
      PARAMETER (FOUR=4.0D+00)
      PARAMETER (PT5 =0.5D+00)
      DIMENSION X(*),W(*)
      PI=FOUR*ATAN(ONE)
      NR=(N+1)/2
      XM=PT5*(X2+X1)
      XL=PT5*(X2-X1)
      PIN=(PI/TWO)/(2*N+1)
      DO I=1,NR
         Z=COS(PIN*(4*I-1))
    1    CONTINUE
            P1=ONE
            P2=ZERO
            DO J=1,N
               P3=P2
               P2=P1
               P1=((J+J-1)*Z*P2-(J-1)*P3)/J
            ENDDO
            PP=N*(Z*P1-P2)/(Z*Z-ONE)
            Z1=Z
            Z =Z1-P1/PP
            IF(ABS(Z-Z1).GT.EPS) GO TO 1
         X(    I)=XM+XL*Z
         X(N+1-I)=XM-XL*Z
         W(    I)=TWO*XL/((ONE-Z*Z)*PP*PP)
         W(N+1-I)=W(I)
      ENDDO
      RETURN
      END
C*MODULE DFTGRD  *DECK GRDDFT
      SUBROUTINE GRDDFT(L2,LAST)
C***********************************************************************
C
C     THIS PROGRAM UTILIZES AN ARBITRARY DENSITY TO CALCULATE
C     THE CORRELATION CORRECTION TO SCF WITH AN ARBITRARY
C     SET OF DENSITY FUNCTIONALS.
C
C     --- WRITTEN BY ---
C     MUNEAKI KAMIYA, TAKAO TSUNEDA, SUSUMU YANAGISAWA, DMITRI FEDOROV
C
C                                                UNIVERSITY OF TOKYO
C
C     THE ORIGINAL PROGRAM WAS WRITTEN BY
C     NEVIN OLIPHANT, QUANTUM THEORY PROJECT, UNIVERSITY OF FLORIDA
C     HIDEO SEKINO, QUANTUM THEORY PROJECT, UNIVERSITY OF FLORIDA
C
C***********************************************************************
C
C     THIS ROUTINE CONSISTS OF FOUR PARTS
C       1. MOMORY ALLOCATION
C       2. GEOMETRY AND SYMMETRY SETTING                        (DFTSET)
C       3. CALCULATING EXCHANGE CORRELATION ENERGY            (DFTEXCOR)
C       4. CALCULATING EXCHANGE CORRELATION ENERGY GRADIENT    (DFTGRAD)
C
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*10 XCHNG,CRRN,PFCHR
      CHARACTER*3 AGROUP
C     ----
      LOGICAL UROHF,OPTGRD,GOPARR,DSKWRK,MASWRK,NEEDGR,SAVEGRID,SVDSKW,
     *        ABELPT,GRDOUT
C     ---- MEMORY ADDRESS -----
      SAVE ISTART,IPCOEFF,IEXPS,IWGHT,IXDAT,IYDAT,IZDAT,
     *     IATMXVC,IATMYVC,IATMZVC,IRI,IRIJ,IWTINTR,IAOX,IGAOX,IGAOY,
     *     IGAOZ,ITOTWT,IRSQRD,IAIJ,IGLROOT,IGLWGHT,IANGXV,IANGYV,IANGZV
     *    ,IPTRAD,IWTRAD,IVMOA,IDVMOXA,IDVMOYA,IDVMOZA,IVALGA,IVMOB,
     *     IDVMOXB,IDVMOYB,IDVMOZB,IVALGB,IFACTR,INATM,IFACT,NAPTS,
     *     IANGN,INEQATM,INDEG,IUNIQUE
C
      SAVE ISYMXY,ISYMXZ,ISYMYZ,ISYMRX,ISYMRY,ISYMRZ,ISYMI,
     *     ISYMC1,NANGPT,NSYMAT,NANG,NDER,NEEDGR,NDFTEN
C
      SAVE UROHF
C     -------
      PARAMETER (MXATM=2000,MXGTOT=20000)
      PARAMETER (ZERO=0.0D+00,BOHR=0.52917715D+00)
C
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      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 /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DIMENSION X(*),FA(*),FB(*),COEFFA(*),COEFFB(*),G2AXX(*),G2AYY(*),
     *          G2AZZ(*),G2AXY(*),G2AYZ(*),G2AXZ(*),DA(*),DB(*),DEDFT(*)
     *         ,UVEC(*),DWTINT(*),DWTTOT(*),DTOTWT(*)
      DIMENSION IJX(35),IJY(35),IJZ(35),BSLRD(137)
C
      DATA UHF,ROHF/8HUHF     ,8HROHF    /
C***********************************************************************
C     BRAGG-SLATER RADII FOR DETERMINING THE RELATIVE SIZE OF THE
C     POLYHEDRA IN THE POLYATOMIC INTEGRATION SCHEME
C***********************************************************************
      DATA BSLRD  /0.529177D+00,0.31D+00,1.45D+00,1.05D+00,0.85D+00,
     1             0.70D+00, 0.65D+00, 0.60D+00, 0.50D+00, 0.38D+00,
     2             1.80D+00, 1.50D+00, 1.25D+00, 1.10D+00, 1.00D+00,
     3             1.00D+00, 1.00D+00, 0.71D+00, 2.20D+00, 1.80D+00,
     4             1.60D+00, 1.40D+00, 1.35D+00, 1.40D+00, 1.40D+00,
     5             1.40D+00, 1.35D+00, 1.35D+00, 1.35D+00, 1.35D+00,
     6             1.30D+00, 1.25D+00, 1.15D+00, 1.15D+00, 1.15D+00,
     7             0.88D+00, 2.35D+00, 2.00D+00, 1.80D+00, 1.55D+00,
     8             1.45D+00, 1.45D+00, 1.35D+00, 1.30D+00, 1.35D+00,
     9             1.40D+00, 1.60D+00, 1.55D+00, 1.55D+00, 1.45D+00,
     1             1.45D+00, 1.40D+00, 1.40D+00, 1.08D+00, 2.60D+00,
     2             2.15D+00, 1.95D+00, 1.85D+00, 1.85D+00, 1.85D+00,
     3             1.85D+00, 1.85D+00, 1.85D+00, 1.80D+00, 1.75D+00,
     4             1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00,
     5             1.75D+00, 1.55D+00, 1.45D+00, 1.35D+00, 1.35D+00,
     6             1.30D+00, 1.35D+00, 1.35D+00, 1.35D+00, 1.50D+00,
     7             1.90D+00, 1.80D+00, 1.60D+00, 1.90D+00, 1.27D+00,
     8             1.20D+00, 2.60D+00, 2.15D+00, 1.95D+00, 1.80D+00,
     9             1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00,
     1             1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00,
     2             1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00,
     2             1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00,
     2             1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00,
     2             1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00,
     2             1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00,
     2             1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00,
     2             1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00, 1.75D+00,
     2             1.75D+00, 1.75D+00/
      DATA IJX    / 1, 2, 1, 1, 3, 1, 1, 2, 2, 1,
     1              4, 1, 1, 3, 3, 2, 1, 2, 1, 2,
     2              5, 1, 1, 4, 4, 2, 1, 2, 1, 3,
     3              3, 1, 3, 2, 2/
      DATA IJY    / 1, 1, 2, 1, 1, 3, 1, 2, 1, 2,
     1              1, 4, 1, 2, 1, 3, 3, 1, 2, 2,
     2              1, 5, 1, 2, 1, 4, 4, 1, 2, 3,
     3              1, 3, 2, 3, 2/
      DATA IJZ    / 1, 1, 1, 2, 1, 1, 3, 1, 2, 2,
     1              1, 1, 4, 1, 2, 1, 2, 3, 3, 2,
     2              1, 1, 5, 1, 2, 1, 2, 4, 4, 1,
     3              3, 3, 2, 2, 3/
C
      UROHF = SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
C
C     FOR SLATER EXHANGE AND LDA OR NO CORRELATION DO NO NEET TO HANDLE
C     THE GRADIENT TERMS
C
      NEEDGR=.NOT.(NEXFG.EQ.1.AND.(NCORFG.EQ.0.OR.NCORFG.EQ.1))
      PI=ATAN(1.0D+00)*4.0D+00
C
C*********************************************************************
C     ----- GET GRID DFT MEMORY -----
C*********************************************************************
C     IN THE FOLLOWING, ARRAYS ARE ALLOCATED AS NAT, WHEREAS SOME ARE
C     IN FACT USED AS NSYMAT. THIS MAY SAVE MEMORY BUT REQUIRES SMALL
C     CHANGES IN THE CODE.
C
C     NXT = IBTYP.EQ.1
      NANGPT = NTHE*NPHI
      NPT = NRAD*NANGPT
      NPTME=(NPT-1)/NPROC+1
      IF(IBTYP.EQ.1) NPTME=NPT
C     TOTWT CANNOT BE DIVIDED FOR DYNAMIC LOAD BALANCING SINCE WEIGHTS
C     WILL BE REQUIRED AT "RANDOM" (DDI ARRAYS MAY WORK).
      CALL BASCHK(MAXL)
      CALL DERCHK(NDER)
      NANG=MAXL+1+1
C     1 FOR GRADIENT TERMS
      NANG1=NANG+1+NDER
C
C     ----- FOR SUBROUTINE INPBAS -----
C
      IPCOEFF = 1
      IEXPS   = IPCOEFF + MXGTOT
      ISTART  = IEXPS   + MXGTOT
C
      INDEG   = ISTART
      INEQATM = INDEG   + (NAT-1)/NWDVAR+1
      IUNIQUE = INEQATM + (NAT*NAT-1)/NWDVAR+1
      ISTART  = IUNIQUE + (NAT-1)/NWDVAR+1
C
C     ----- FOR ROUTINE GRDPT -----
C
      IWGHT   = ISTART
      IXDAT   = IWGHT   + NANGPT*NAT
      IYDAT   = IXDAT   + NANGPT*NAT
      IZDAT   = IYDAT   + NANGPT*NAT
      IATMXVC = IZDAT   + NANGPT*NAT
      IATMYVC = IATMXVC + NAT*NAT
      IATMZVC = IATMYVC + NAT*NAT
      IRI     = IATMZVC + NAT*NAT
      IRIJ    = IRI     + NAT
      IWTINTR = IRIJ    + NAT*NAT
      IRSQRD  = IWTINTR + NAT
      IAIJ    = IRSQRD  + NAT
      IGLROOT = IAIJ    + NAT*NAT
      IGLWGHT = IGLROOT + NTHE*NTHE
      IPTRAD  = IGLWGHT + NTHE*NTHE
      IWTRAD  = IPTRAD  + NRAD
      IANGXV  = IWTRAD  + NRAD
      IANGYV  = IANGXV  + NAT*NANG1
      IANGZV  = IANGYV  + NAT*NANG1
      IFACTR  = IANGZV  + NAT*NANG1
      INATM   = IFACTR  + (NAT-1)/NWDVAR+1
      IFACT   = INATM   + (NAT-1)/NWDVAR+1
      NAPTS   = IFACT   + (NAT-1)/NWDVAR+1
      IANGN   = NAPTS   + (NAT-1)/NWDVAR+1
      LAST    = IANGN   + (NAT*2-1)/NWDVAR+1
C
C     ----- FOR SUBROUTINE DMATD ------
C
      IVMOA  = LAST
      IDVMOXA= IVMOA  + NUM
      IDVMOYA= IDVMOXA+ NUM
      IDVMOZA= IDVMOYA+ NUM
      IVMOB  = IDVMOZA+ NUM
      IDVMOXB= IVMOB  + NUM
      IDVMOYB= IDVMOXB+ NUM
      IDVMOZB= IDVMOYB+ NUM
      IVALGA = IDVMOZB+ NUM
      IVALGB = IVALGA + NUM
      LAST   = IVALGB + NUM
C
      ITOTWT  = LAST
      IAOX    = ITOTWT  + NAT*NPTME
      IGAOX   = IAOX    + NUM
      IGAOY   = IGAOX   + NUM
      IGAOZ   = IGAOY   + NUM
      LAST    = IGAOZ   + NUM
C     WRITE(6,*) 'DFTALLOC',NDFTEN,NANGPT,NAT,NTHE,NRAD,NPTME
      NDFTEN=LAST
      RETURN
C***********************************************************************
      ENTRY DFTSET(X,NPRT,OPTGRD)
C***********************************************************************
      NFTDFT=22
      IF(OPTGRD) THEN
C       ISYMC1 MUST HAVE BEEN SET DURING A PRECEDING ENERGY CALCULATION
C       FTNCHECK WILL GIVE A WARNING
        SAVEGRID=ISYMC1.EQ.1.AND.NDER.GT.0
        IF(SAVEGRID) THEN
          SVDSKW=DSKWRK
          DSKWRK=.TRUE.
          CALL SEQOPN(NFTDFT,'DFTGRID','UNKNOWN',.FALSE.,'UNFORMATTED')
          READ(NFTDFT) (X(I),I=1,NDFTEN)
          CALL SEQCLO(NFTDFT,'KEEP')
          DSKWRK=SVDSKW
          RETURN
C         READ THE GRID AND RETURN
        ENDIF
      ENDIF
C     ---------------------------------
C           WRITE FUNCTIONAL NAME IF NPRT.NE.0
C     ---------------------------------
C
      NPT=NRAD*NTHE*NPHI
      IF(NPRT.NE.0.AND.MASWRK) THEN
        CALL RNAMEXC(XCHNG,CRRN)
        WRITE(IW,1010) XCHNG
        WRITE(IW,1020) CRRN
        IF (NEXFG.EQ.6) THEN
          CALL REFPFREE (PFCHR)
          WRITE(IW,1030) PFCHR
        ENDIF
        DFTTHR0=1.0D-03/(NPT*NAT)
        IF(DFTTHR.NE.ZERO) DFTTHR0=DFTTHR
        IF(DFTTHR0.LT.1.1D-15) THEN
          WRITE(IW,1060)
        ELSE
          WRITE(IW,1040) DFTTHR0
        ENDIF
        IF(NRAD*NPHI*NTHE.NE.NRAD0*NPHI0*NTHE0) WRITE(IW,1050) SW0
      ENDIF
 1010 FORMAT(/5X,'EXCHANGE FUNCTIONAL   =',A10)
 1020 FORMAT(5X, 'CORRELATION FUNCTIONAL=',A10)
 1030 FORMAT(5X, 'KINETIC FUNCTIONAL    =',A10)
 1040 FORMAT(5X, 'DFT THRESHOLD         =',E8.3)
 1050 FORMAT(5X, 'GRID CHANGE THRESHOLD =',E8.3)
 1060 FORMAT(5X, 'ALL DFT THRESHOLDS ARE TURNED OFF.')
 1070 FORMAT(5X, 'USING ',A3,' ABELIAN SUBGROUP FOR OCTANTS AND',
     *           ' FULL SYMMETRY FOR ATOMS.',/)
C
C     ---------------------------------------------------------
C           GET SYMMETRY AND INITIAL SYMMETRY INFORMATION
C     ---------------------------------------------------------
C
      CALL INPINF(ISYMXY,ISYMXZ,ISYMYZ,ISYMRX,ISYMRY,ISYMRZ,ISYMI,
     >            ISYMC1,AGROUP)
      IF(NPRT.NE.0.AND.AGROUP.NE.'C1'.AND..NOT.ABELPT().AND.MASWRK)
     *   WRITE(IW,1070) AGROUP
C
C     THERE IS A SUBTLE DIFFERENCE BETWEEN RUNNING GRADIENT AND LOCAL
C     FRAME, BOTH WITHOUT OCTANT SYMMETRY.
C     FOR THE GRADIENT, WE GENERATE GRID WITH THE OCTANT SYMMETRY AND
C     USE IT FOR THE ENERGY, BUT NOT FOR THE GRADIENT.
C     FOR THE LOCAL FRAME RUNS, DO NOT USE THE OCTANT SYMMETRY.
C     THIS IS WHY LOCFRAME IS PROVIDED AS AN ARGUMENT TO INPINF.
C
C     ------------------------------------------------
C           SET UP PRIMITIVE COEFFICIENT VECTORS
C     ------------------------------------------------
C
      CALL SETPNRM
      CALL INPBAS(X(IPCOEFF))
C
C     --------------------------------------------------------------
C                 THE SETUP OF THE FACTOR IFACTR
C           FOR THE SYMMETRY-UNIQUE ATOM TO BE MULTIPLIED BY.
C     --------------------------------------------------------------
C
      CALL SYMUNQ(X(IFACTR),X(INATM),NSYMAT,X(INDEG),X(INEQATM),
     *            X(IUNIQUE),IZCORE)
C
      CALL ATMVEC(X(IATMXVC),X(IATMYVC),X(IATMZVC),X(IRIJ))
C
      CALL RADPT(X(IPTRAD),X(IWTRAD),NRAD)
      IF(IBTYP.EQ.1) CALL VCLR(X(ITOTWT),1,NAT*NPT)
C
      DO 30 NCNTR = 1, NAT
C
C     ---- DUMMY, GHOST AND SYMMETRY NON-UNIQUE ATOMS ARE SKEPT. ----
C
        IF (IXFTCH(X(IUNIQUE),NCNTR).EQ.0) GOTO 30
C
C     ------------------------------------------------
C           DETERMINE WHICH QUADRANTS ARE UNIQUE
C                 AND THE SYMMETRY FACTOR
C              TO MULTIPLY FINAL INTEGRAL BY.
C     -------------------------------------------------
C
        CALL SYMFCT(ISYMYZ,ISYMXZ,ISYMXY,ISYMRX,ISYMRY,
     >              ISYMRZ,ISYMI,NCNTR,
     >              IQ1,IQ2,IQ3,IQ4,IQ5,IQ6,IQ7,IQ8,
     >              X(IFACT),ISYMC1)
C
C     ----------------------------------------------------
C            DETERMINE THE GRID POINTS FOR NCNTR
C                THE FUZZY CELL METHOD OF BECKE
C           ( A.D.BECKE, J.CHEM.PHYS.,88,2547,1988).
C     ----------------------------------------------------
C
C
        NCNX = NCNTR
        CALL GRDPT(NCNX,X(IWGHT),X(IXDAT),X(IYDAT),X(IZDAT),X(IATMXVC),
     *             X(IATMYVC),X(IATMZVC),X(IRI),X(IRIJ),X(IWTINTR),
     *             X(ITOTWT),X(IAIJ),X(IGLROOT),X(IGLWGHT),BSLRD,
     *             X(IANGN),X(IFACT),X(NAPTS),X(IPTRAD),X(IWTRAD),IQ2,
     *             IQ3,IQ4,IQ5,IQ6,IQ7,NRAD,NTHE,NPHI,OPTGRD,
     *             BOHR,IZCORE)
 30   CONTINUE
      IF(IBTYP.EQ.1) CALL DDI_GSUMF(2314,X(ITOTWT),NAT*NPT)
      SAVEGRID=ISYMC1.EQ.1.AND.NDER.GT.0.AND.
     *         NRAD*NPHI*NTHE.GE.NRAD0*NPHI0*NTHE0
C     THE LATTER CONDITION SAVES SOME TIME BY NOT WRITING THE COARSER
C     GRID TO DISK SINCE WE SHALL LATER DEFINE A FINER GRID.
      IF(SAVEGRID) THEN
C        SAVE THE GRID FOR THE GRADIENT AND HOPEFULLY THE HESSIAN
         SVDSKW=DSKWRK
         DSKWRK=.TRUE.
         CALL SEQOPN(NFTDFT,'DFTGRID','UNKNOWN',.FALSE.,'UNFORMATTED')
         WRITE(NFTDFT) (X(I),I=1,NDFTEN)
         CALL SEQCLO(NFTDFT,'KEEP')
         DSKWRK=SVDSKW
C        CALL DAWRIT(IDAF,IODA,X,NDFTEN,400,0)
      ENDIF
      RETURN
C***********************************************************************
      ENTRY DFTEXCOR(X,FA,FB,COEFFA,COEFFB,L1,L2,EEXC,TOTELE)
C***********************************************************************
C
C     COMPUTE THE EXCHANGE-CORRELATION ENERGY FUNCTIONAL PART
C
      CALL VCLR(FA,1,L2)
      IF (UROHF) THEN
         CALL VCLR(FB,1,L2)
      ENDIF
C
C     ----- LOOP FOR ATOMS -----
C
      TOTELE = ZERO
      EEXC   = ZERO
      DO 930 NCNTR = 1, NAT
        INC0=IXFTCH(X(IUNIQUE),NCNTR)
        IF(INC0.EQ.0) GOTO 930
C
C***********************************************************************
C
C     ----------------------------------------------------------------
C            CALCULATE THE EXCHANGE-CORRELATION ENERGY FUNCTIONAL
C     ----------------------------------------------------------------
C
        CALL DMATD(X(ITOTWT),X(IFACT),X(NAPTS),COEFFA,COEFFB,X(IANGN),
     *             X(IFACTR),INC0,X(IVMOA),X(IDVMOXA),X(IDVMOYA),
     *             X(IDVMOZA),X(IVALGA),X(IVMOB),X(IDVMOXB),X(IDVMOYB),
     *             X(IDVMOZB),X(IVALGB),FA,FB,NCNTR,EEXC,TOTELE,NANG,
     *             X(IPTRAD),X(IXDAT),X(IYDAT),X(IZDAT),X(IATMXVC),
     *             X(IATMYVC),X(IATMZVC),X(IANGXV),X(IANGYV),X(IANGZV),
     *             BSLRD,X(IRSQRD),X(IPCOEFF),X(IEXPS),X(IAOX),X(IGAOX),
     *             X(IGAOY),X(IGAOZ),IJX,IJY,IJZ,L1,BOHR,NEEDGR,IZCORE,
     *             UROHF)
C
 930  CONTINUE
C     CALL DDI_GSUMF(2310,FA,L2)
C     IF(UROHF) CALL DDI_GSUMF(2311,FB,L2)
C     CALL DDI_GSUMF(2312,EEXC,1)
C     CALL DDI_GSUMF(2313,TOTELE,1)
      RETURN
C***********************************************************************
      ENTRY DFTGRAD(X,G2AXX,G2AYY,G2AZZ,G2AXY,G2AYZ,G2AXZ,
     *              COEFFA,COEFFB,DA,DB,DEDFT,UVEC,DWTINT,DWTTOT,DTOTWT,
     *              L1,TOTELE,GRDOUT)
C***********************************************************************
C
C     ----- LOOP FOR ATOMS -----
C
      TOTELE = ZERO
      DO 940 NCNTR = 1, NAT
        INC0=IXFTCH(X(IUNIQUE),NCNTR)
        IF (INC0.EQ.0) GOTO 940
        NCNX = NCNTR
        CALL DEREXC(X(ITOTWT),X(NAPTS),COEFFA,COEFFB,DA,DB,X(IANGN),
     *              X(IFACTR),INC0,X(IVMOA),X(IDVMOXA),X(IDVMOYA),
     *              X(IDVMOZA),X(IVMOB),X(IDVMOXB),X(IDVMOYB),X(IDVMOZB)
     *             ,X(IWGHT),X(IRI),X(IRIJ),X(IWTINTR),X(IAIJ),X(IWTRAD)
     *             ,NCNX,L1,TOTELE,NANG+NDER,X(IPTRAD),X(IXDAT),
     *              X(IYDAT),X(IZDAT),X(IATMXVC),X(IATMYVC),X(IATMZVC),
     *              X(IANGXV),X(IANGYV),X(IANGZV),BSLRD,X(IRSQRD),
     *              X(IPCOEFF),X(IEXPS),X(IAOX),X(IGAOX),X(IGAOY),
     *              X(IGAOZ),G2AXX,G2AYY,G2AZZ,G2AXY,G2AYZ,G2AXZ,IJX,
     *              IJY,IJZ,BOHR,NEEDGR,DEDFT,IZCORE,UVEC,DWTINT,DWTTOT,
     *              DTOTWT,GRDOUT)
 940  CONTINUE
      RETURN
C1193 FORMAT(/'TOTAL EXC ENERGY',F20.10,'TOTAL ELECTRON',F20.10)
      END
C*MODULE DFTGRD  *DECK GRDOCT
      SUBROUTINE GRDOCT(NAT,NITR,NCNTR,ATMXVEC,ATMYVEC,ATMZVEC,RI,RIJ,
     *                  AIJ,WGHT,RADWT,P,UVEC,DP,DZ,DW)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00,TWO=2.0D+00)
C
C     COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,NDFTFG,
C    *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0
C     COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION ATMXVEC(NAT,NAT),ATMYVEC(NAT,NAT),ATMZVEC(NAT,NAT),
     *          RI(NAT),RIJ(NAT,NAT),P(NAT),AIJ(NAT,NAT),
     *          UVEC(3,NAT),DMUJI(3),DP(3,NAT,NAT),DZ(3,NAT),DW(3,NAT)
C**********************************************************************
C
C     COMPUTE WEIGHT(=CELL VOLUME) DERIVATIVES.
C     THE FORMULAE AND NOTATION ARE TAKEN FROM
C     JOHNSON ET AL, JCP 98, 5612 (1993)
C
C     THE RESULS ARE PUT INTO DW.
C**********************************************************************
C
C     THRESHOLDS FOR THE REMARK 62 IN THE ABOVE REFERENCE. IT
C     CORRESPONDS TO A*B/B, WHEN B IS ZERO. DIVISION BY B IS AVOIDED
C     FOR B LESS THAN THRETH (AND THUS ASSIGNING A THE VALUE OF ZERO).
C     THIS RELIES ON A BEING ZERO WHENEVER B IS ZERO, WHICH APPEARS TO
C     BE TRUE ALGEBRAICALLY. IF |A|>THRETH1, THEN THE CODE ABORTS.
C     IF IT IS FOUND THAT THIS HAPPENS, IT WOULD BE NECESSARY TO
C     CALCULATE A EXPLICITLY AND AVOID *B/B, WHICH, HOWEVER, IS A BIT
C     OF WORK.  (NOTE G4 BELOW IS NOT A, BUT A IS PROPORTIONAL TO G4).
C
      THRETH=1.0D-12
      THRETH1=1.0D-08
C
      CALL DACOPY(NAT,ONE,P,1)
      CALL VCLR(DP,1,3*NAT*NAT)
      CALL VCLR(DZ,1,3*NAT)
      CALL VCLR(DW,1,3*NAT)
C
C        LOOP JATM RUNS OVER THE DERIVATIVE COORDINATES B
C        LOOP IATM RUNS OVER CENTRES A
C
         DO 611 JATM=1,NAT
            IF (ABS(AIJ(JATM,1)+1.0D+00).LT.1.0D-05) THEN
              P(JATM) = 0.0D+00
              GOTO 611
            ENDIF
            DO 612 IATM=1,NAT
               IF(IATM.EQ.JATM) GOTO 612
               IF(ABS(AIJ(JATM,IATM)-1.0D+00).LT.1.0D-05) GOTO 612
               RBA=RIJ(JATM,IATM)
C
C     COMPUTE HYPERBOLIC COORDINATES MU(B,A) (XMUJI)
C     AND MU'(B,A)=MU(B,A)+A(B,A)*(1-MU(B,A)**2)
C
               ZMUJI=(RI(JATM)-RI(IATM))/RBA
               XMUJI=ZMUJI+AIJ(JATM,IATM)*(ONE-ZMUJI*ZMUJI)
C              XMUJI=ZMUJI
C
C     BECKE'S FUZZY CELL METHOD FOR MOLECULAR GRID QUADRATURE
C       P1(MU)=P(MU)=3/2*MU-1/2*MU**3, PN(MU)=P(P{N-1}(MU))
C       CUTIJ=S(MUJI)=1/2*(1-PN(MU))
C
               F4=XMUJI
               G4=ONE
C
               DO K=1,NITR
                  G4= G4*(ONE-F4*F4)
                  F4= F4*(1.5D+00-0.5D+00*F4*F4)
               END DO
C
C              COMPUTE S(MU(A,B)) AND S(MU(B,A))
C
               F2   = 0.5D+00*F4
               CUTIJ= 0.5D+00+F2
               CUTJI= 0.5D+00-F2
               P(JATM)=P(JATM)*CUTJI
C
               IF(JATM.EQ.NCNTR) GOTO 612
C
C              THE CODE BELOW DOES NOT WORK FOR GRAD(B) X, FOR SOME X,
C              IN CASE OF B==NCNTR, DUE TO A MORE COMPLEX DERIVATIVE,
C              SO WE USE TRANSLATIONAL INVARIANCE INSTEAD
C              (X IS MU(B,A), P(B), W(A)).
C              THE LINES ASSIGNING ID WILL MAKE GRAD(B) MU(B,A) CORRECT
C              EVEN FOR B==NCNTR (WHICH WE DO NOT NORMALLY NEED).
C              GRAD(B) P(A) WILL HOWEVER STILL BE INCORRECT (B==NCNTR).
C              IF(JATM.EQ.NCNTR) THEN
C                 ID=IATM
C              ELSE
C                 ID=JATM
C              ENDIF
               ID=JATM
C
C              COMPUTE GRAD(B) MU'(B,A) = A * GRAD(B) MU(B,A)
C                      A = 1-2*A(B,A)*MU(B,A)
C
               A=ONE-TWO*AIJ(JATM,IATM)*ZMUJI
               B=ZMUJI/(RBA**2)
               DMUJI(1)=A*(-UVEC(1,ID)/RBA-B*ATMXVEC(JATM,IATM))
               DMUJI(2)=A*(-UVEC(2,ID)/RBA-B*ATMYVEC(JATM,IATM))
               DMUJI(3)=A*(-UVEC(3,ID)/RBA-B*ATMZVEC(JATM,IATM))
C              WRITE(6,*) 'MU=BA',XMUJI,'GRAD(B) MU(B,A)=',JATM,IATM,
C    >                    DMUJI(1),DMUJI(2),DMUJI(3)
C
C              COMPUTE D S(MU(B,A))/D MU(B,A)
C              THE COEFFICIENT IN DCOEFF=-1/2*(3/2)**NITR
C
               DCOEFF = -81.0D+00/32.0D+00*G4
C              DCOEFF = -27.0D+00/16.0D+00*G4
C
C              COMPUTE 1/P(A)*GRAD(B) P(A) (DP(1,IATM,JATM))
C
               IF (ABS(CUTIJ).GT.THRETH) THEN
                  DCUTIJ = DCOEFF/CUTIJ
                  DP(1,IATM,JATM)= -DCUTIJ*DMUJI(1)
                  DP(2,IATM,JATM)= -DCUTIJ*DMUJI(2)
                  DP(3,IATM,JATM)= -DCUTIJ*DMUJI(3)
              ELSE
                 IF(ABS(G4).GT.THRETH1) THEN
                    WRITE(6,*) 'G4',G4
                    CALL ABRT
                 ENDIF
              ENDIF
              IF (ABS(CUTJI).GT.THRETH) THEN
                 DCUTJI = DCOEFF/CUTJI
                 DP(1,JATM,JATM)=DP(1,JATM,JATM)+DCUTJI*DMUJI(1)
                 DP(2,JATM,JATM)=DP(2,JATM,JATM)+DCUTJI*DMUJI(2)
                 DP(3,JATM,JATM)=DP(3,JATM,JATM)+DCUTJI*DMUJI(3)
              ELSE
                 IF(ABS(G4).GT.THRETH1) THEN
                    WRITE(6,*) 'G4',G4
                    CALL ABRT
                 ENDIF
              ENDIF
 612       CONTINUE
 611    CONTINUE
C
C       COMPUTE Z (TOTAL WEIGHT)
C
        Z=ZERO
        DO IATM=1,NAT
          Z=Z+P(IATM)
        ENDDO
C       WRITE(6,*) '=Z',Z
C       ATMWT=P(NCNTR)/Z
        SPHWT=RADWT*WGHT
C       TOTWT=ATMWT*SPHWT
C       WRITE(6,*) 'W(A)',NCNTR,TOTWT
C       DO IATM=1,NAT
C         WRITE(6,*) 'P(A)',IATM,P(IATM)
C          DO JATM=1,NAT
C           WRITE(6,*) NCNTR,'=GRAD(B) P(A)',JATM,IATM,DP(1,IATM,JATM)
C    *        *P(IATM),DP(2,IATM,JATM)*P(IATM),DP(3,IATM,JATM)*P(IATM)
C         ENDDO
C       ENDDO
        WTC=P(NCNTR)/Z*SPHWT
C       IF(ABS(WTC).GT.WDCUTOFF) THEN
        DO 690 JATM=1,NAT
           IF (JATM.EQ.NCNTR) GOTO 690
C
C          COMPUTE GRAD(B) Z (DZ(*,JATM))
C
           DO IATM=1,NAT
              DZ(1,JATM)=DZ(1,JATM)+DP(1,IATM,JATM)*P(IATM)
              DZ(2,JATM)=DZ(2,JATM)+DP(2,IATM,JATM)*P(IATM)
              DZ(3,JATM)=DZ(3,JATM)+DP(3,IATM,JATM)*P(IATM)
           ENDDO
C          WRITE(6,*) '=GRAD(B) Z',JATM,DZ(1,JATM),DZ(2,JATM),
C    *                 DZ(3,JATM)
C
C          COMPUTE GRAD(B) W(NCNTR) (DW(*,JATM))
C          DERIVATIVE OF CELL VOLUMES
C
           DW(1,JATM)=WTC*(DP(1,NCNTR,JATM)-DZ(1,JATM)/Z)
           DW(2,JATM)=WTC*(DP(2,NCNTR,JATM)-DZ(2,JATM)/Z)
           DW(3,JATM)=WTC*(DP(3,NCNTR,JATM)-DZ(3,JATM)/Z)
C          WRITE(6,*) '==GRAD(B) W(A)',JATM,NCNTR,DW(1,JATM),
C    *                 DW(2,JATM),DW(3,JATM)
           DW(1,NCNTR)=DW(1,NCNTR)-DW(1,JATM)
           DW(2,NCNTR)=DW(2,NCNTR)-DW(2,JATM)
           DW(3,NCNTR)=DW(3,NCNTR)-DW(3,JATM)
 690    CONTINUE
C       ENDIF
C       DO JATM=1,NAT
C          WRITE(6,*) '==GRAD(B) W(A)',JATM,NCNTR,DW(1,JATM),
C    *                 DW(2,JATM),DW(3,JATM)
C       ENDDO
C
      RETURN
      END
C*MODULE DFTGRD  *DECK GRDPT
      SUBROUTINE GRDPT(NCNTR,WGHT,XDAT,YDAT,ZDAT,ATMXVEC,ATMYVEC,ATMZVEC
     *                ,RI,RIJ,WTINTR,TOTWT,AIJ,GLROOT,GLWGHT,BSLRDII,
     *                 IANGN,IIFACT,NAPTS,PTRAD,WTRAD,IQ2,IQ3,IQ4,IQ5,
     *                 IQ6,IQ7,NRADPT,NTHE,NPHI,OPTGRD,BOHR,IZCORE)
C***********************************************************************
C
C     THIS ROUTINE SETS UP THE NUMERICAL GRID FOR DIFFERENT
C     ABELIAN SYMMETRIES.
C
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL   OPTGRD,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000,ONE=1.0D+00,TWO=2.0D+00,HALF=0.5D+00)
C
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /DFTEXC/ PI,CSLT,CB88,CLYP,CVWN,QOP,NEXFG,NCORFG,NPFFG
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION WGHT(NAT,NTHE*NPHI),XDAT(NAT,NTHE*NPHI),
     >          YDAT(NAT,NTHE*NPHI),ZDAT(NAT,NTHE*NPHI),
     >          ATMXVEC(NAT,NAT),ATMYVEC(NAT,NAT),ATMZVEC(NAT,NAT),
     >          RI(NAT),RIJ(NAT,NAT),WTINTR(NAT),TOTWT(NAT,*),
     >          AIJ(NAT,NAT),GLROOT(NTHE,NTHE),GLWGHT(NTHE,NTHE),
     *          BSLRDII(137),IANGN(NAT,2),IIFACT(NAT),NAPTS(NAT),
     *          PTRAD(NRADPT),WTRAD(NRADPT),IZCORE(*)
C
C***********************************************************************
C
C     EULER-MACLAURIN SCHEME, RADIAL INTEGRATION
C     GAUSS-MARKOV QUADRATURE, ANGULAR INTEGRATION
C     BECKE POLYATOMIC SCHEME
C
C***********************************************************************
C     THE # OF ITERATIONS OF THE EQUATION WHICH CREATES THE 'FUZZY'
C     BOUNDARY
C***********************************************************************
      NITR = 4
C***********************************************************************
C
C     SET UP ATOMIC VECTORS AND CALCULATE THE SURFACE SHIFTING
C     PARAMETER AIJ
C
C***********************************************************************
      IFACT = IIFACT(NCNTR)
      DO 6010 IATM=1,NAT
         AIJ(IATM,IATM) = 0.0D+00
         IBSR1=INT(ZAN(IATM)+IZCORE(IATM))
         IF(NFG.NE.0) IBSR1=IAN(IATM)
         IF(IBSR1.EQ.0) THEN
           AIJ(IATM,1)=-1.0D+00
           GOTO 6010
         ENDIF
         DO 6020 JATM=1,IATM-1
            IBSR2=INT(ZAN(JATM)+IZCORE(JATM))
            IF(NFG.NE.0) IBSR2=IAN(JATM)
            IF(IBSR2.EQ.0) THEN
              AIJ(IATM,JATM)= 1.0D+00
              GOTO 6020
            ENDIF
            CHI=BSLRDII(IBSR1)/BSLRDII(IBSR2)
            CHI2=(CHI-ONE)/(CHI+ONE)
            AIJ(IATM,JATM)=CHI2/(CHI2*CHI2-ONE)
            IF(AIJ(IATM,JATM).GT.0.5D+00) AIJ(IATM,JATM)=0.5D+00
            IF(AIJ(IATM,JATM).LT.-0.5D+00) AIJ(IATM,JATM)=-0.5D+00
 6020 CONTINUE
         DO 6030 JATM=IATM+1,NAT
            IBSR2=INT(ZAN(JATM)+IZCORE(JATM))
            IF(NFG.NE.0) IBSR2=IAN(JATM)
            IF(IBSR2.EQ.0) THEN
              AIJ(IATM,JATM)= 1.0D+00
              GOTO 6030
            ENDIF
            CHI=BSLRDII(IBSR1)/BSLRDII(IBSR2)
            CHI2=(CHI-ONE)/(CHI+ONE)
            AIJ(IATM,JATM)=CHI2/(CHI2*CHI2-ONE)
            IF(AIJ(IATM,JATM).GT.0.5D+00) AIJ(IATM,JATM)=0.5D+00
            IF(AIJ(IATM,JATM).LT.-0.5D+00) AIJ(IATM,JATM)=-0.5D+00
 6030 CONTINUE
 6010 CONTINUE
C***********************************************************************
C
C     NUMBER OF RADIAL POINTS
C
C***********************************************************************
      ITMP=NRADPT
C***********************************************************************
C
C     SET RADIAL GRID SCALING USING SLATERS RULES
C
C***********************************************************************
      INUC=INT(ZAN(NCNTR)+IZCORE(NCNTR))
      IF(NFG.NE.0) INUC=IAN(NCNTR)
      RAD=BSLRDII(INUC)/BOHR
C
C***********************************************************************
C
C     WEIGHTS FOR THE ANGULAR QUADRATURE ON PRUNE GRID
C
C     GAUSS-LEGENDRE ROOTS AND WEIGHTS FOR THETA GRID
C
C***********************************************************************
      DO ITHE=1,NTHE
      CALL GLGQUD(-ONE,+ONE,GLROOT(1,ITHE),GLWGHT(1,ITHE),ITHE)
      ENDDO
      NTHE=MAX(NTHE,2)
      IF(OPTGRD) IFACT=1
      IF(IFACT.EQ.8)THEN
        ITHEA=1
        ITHEB=NTHE/2
        IPHIA=1
        IPHIB=NPHI/4
        NAPTS(NCNTR)=(NTHE*NPHI)/8
        LTHE=2
        LPHI=4
      ELSE IF(IFACT.EQ.4)THEN
        IF(IQ2.EQ.1)THEN
          ITHEA=1
          ITHEB=NTHE
          IPHIA=1
          IPHIB=NPHI/4
          NAPTS(NCNTR)=(NTHE*NPHI)/4
          LTHE=1
          LPHI=4
        ELSE IF(IQ3.EQ.1)THEN
          ITHEA=1
          ITHEB=NTHE/2
          IPHIA=1
          IPHIA2=3*NPHI/4+1
          IPHIB=NPHI/4
          IPHIB2=NPHI
          NAPTS(NCNTR)=(NTHE*NPHI)/4
          LTHE=2
          LPHI=4
        ELSE IF(IQ5.EQ.1)THEN
          ITHEA=1
          ITHEB=NTHE/2
          IPHIA=1
          IPHIB=NPHI/2
          NAPTS(NCNTR)=(NTHE*NPHI)/4
          LTHE=2
          LPHI=2
        ENDIF
      ELSE IF(IFACT.EQ.2)THEN
        IF(IQ4.EQ.1)THEN
          ITHEA=1
          ITHEB=NTHE
          IPHIA=1
          IPHIA2=3*NPHI/4+1
          IPHIB=NPHI/4
          IPHIB2=NPHI
          NAPTS(NCNTR)=(NTHE*NPHI)/2
          LTHE=1
          LPHI=4
        ELSE IF(IQ6.EQ.1)THEN
          ITHEA=1
          ITHEB=NTHE
          IPHIA=1
          IPHIB=NPHI/2
          NAPTS(NCNTR)=(NTHE*NPHI)/2
          LTHE=1
          LPHI=2
        ELSE IF(IQ7.EQ.1)THEN
          ITHEA=1
          ITHEB=NTHE/2
          IPHIA=1
          IPHIB=NPHI
          NAPTS(NCNTR)=(NTHE*NPHI)/2
          LTHE=2
          LPHI=1
        ENDIF
      ELSE
C       IFACT=1
        ITHEA=1
        ITHEB=NTHE
        IPHIA=1
        IPHIB=NPHI
        NAPTS(NCNTR)=NTHE*NPHI
        LTHE=1
        LPHI=1
      ENDIF
      IF(MOD(NTHE,LTHE).NE.0.OR.MOD(NPHI,LPHI).NE.0) THEN
        IF(MASWRK) WRITE(IW,1000) LTHE,LPHI
C     THIS IS A DOUBLE CHECK.
C     NTHE AND NPHI MUST HAVE ALREADY BEEN ADJUSTED
        CALL ABRT
      ENDIF
      DO ITHE=ITHEA,ITHEB
        WTHE=GLWGHT(ITHE,NTHE)
        COST=GLROOT(ITHE,NTHE)
        SINT=SQRT(ONE-COST**2)
        DO IPHI=IPHIA,IPHIB
          I=(ITHE-ITHEA)*(IPHIB-IPHIA+1)+IPHI
          DPHI=TWO*PI/NPHI
          PHII=DPHI*(IPHI-HALF)
          WPHI=DPHI
          COSP=COS(PHII)
          SINP=SIN(PHII)
          XDAT(NCNTR,I) =SINT*COSP
          YDAT(NCNTR,I) =SINT*SINP
          ZDAT(NCNTR,I) =COST
          WGHT(NCNTR,I) =WTHE*WPHI
        ENDDO
      ENDDO
      LASTI = I
C**********************************************************************
C
C     INTEGRATE OVER ONE OCTANT OF ANGULAR POINTS
C
C**********************************************************************
      IANGN(NCNTR,1)=IPHIA
      IANGN(NCNTR,2)=(ITHEB-ITHEA)*(IPHIB-IPHIA+1)+IPHIB
C
      CALL OCT(ITMP,XDAT,YDAT,ZDAT,ATMXVEC,ATMYVEC,ATMZVEC,RI,RIJ,NAT,
     *         AIJ,WTINTR,TOTWT,WGHT,IANGN,NITR,NCNTR,RAD,NAPTS,
     *         PTRAD,WTRAD)
C**********************************************************************
C     EXTRA DO LOOP WITH REGARD TO IPHI
C                         IN THE CASE OF CI, C2, AND C2H
C**********************************************************************
      IF((IFACT.EQ.4.AND.IQ3.EQ.1).OR.((IFACT.EQ.2.AND.IQ4.EQ.1)))THEN
        GOTO 6050
      ELSE
        RETURN
      ENDIF
 6050 DO ITHE=ITHEA,ITHEB
        WTHE=GLWGHT(ITHE,NTHE)
        COST=GLROOT(ITHE,NTHE)
        SINT=SQRT(ONE-COST**2)
        DO IPHI=IPHIA2,IPHIB2
          I=(ITHE-ITHEA)*(IPHIB2-IPHIA2+1)+IPHI+LASTI-(IPHIA2-1)
          DPHI=TWO*PI/NPHI
          PHII=DPHI*(IPHI-HALF)
          WPHI=DPHI
          COSP=COS(PHII)
          SINP=SIN(PHII)
          XDAT(NCNTR,I) =SINT*COSP
          YDAT(NCNTR,I) =SINT*SINP
          ZDAT(NCNTR,I) =COST
          WGHT(NCNTR,I) =WTHE*WPHI
        ENDDO
      ENDDO
      IANGN(NCNTR,1)=LASTI+1
      IANGN(NCNTR,2)=(ITHEB-ITHEA)*(IPHIB2-IPHIA2+1)+IPHIB2+LASTI
     >    -(IPHIA2-1)
C
      CALL OCT(ITMP,XDAT,YDAT,ZDAT,ATMXVEC,ATMYVEC,ATMZVEC,RI,RIJ,NAT,
     *         AIJ,WTINTR,TOTWT,WGHT,IANGN,NITR,NCNTR,RAD,NAPTS,
     *         PTRAD,WTRAD)
      IANGN(NCNTR,1)=1
      RETURN
 1000 FORMAT(1X,'NTHE HAS TO BE A MULTIPLE OF ',I1,' AND NPHI OF ',I1/,
     *       1X,'IF YOU MUST USE THIS GRID, SET $CONTRL NOSYM=1',/)
      END
C*MODULE DFTGRD  *DECK INPBAS
      SUBROUTINE INPBAS(PCOEFF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXSH=5000, MXGTOT=20000)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      DIMENSION PCOEFF(*)
C
C     FORM A SINGLE PRIMITIVE COEFFICIENT ARRAY
C     USE S FOR L-SHELLS AND TREAT P SEPARATELY
C
      DO I = 1, NSHELL
        DO IK = KSTART(I),KSTART(I)+KNG(I)-1
          IF(KTYPE(I).EQ.1) PCOEFF(IK)=CS(IK)
          IF(KTYPE(I).EQ.2) PCOEFF(IK)=CP(IK)
          IF(KTYPE(I).EQ.2.AND.KMIN(I).EQ.1) PCOEFF(IK)=CS(IK)
          IF(KTYPE(I).EQ.3) PCOEFF(IK)=CD(IK)
          IF(KTYPE(I).EQ.4) PCOEFF(IK)=CF(IK)
          IF(KTYPE(I).EQ.5) PCOEFF(IK)=CG(IK)
        ENDDO
      ENDDO
C
      RETURN
      END
C*MODULE DFTGRD  *DECK INPINF
      SUBROUTINE INPINF(ISYMXY,ISYMXZ,ISYMYZ,ISYMRX,ISYMRY,
     >                  ISYMRZ,ISYMI,ISYMC1,AGROUP)
C***********************************************************************
C
C     THIS SUBROUTINE GETS SYMMETRY AND INITIAL SYMMETRY INFORMATION
C     FOR THE COMPUTATIONAL POINT GROUP.
C
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXATM=2000, ONE=1.0D+00, TM6=1.0D-06)
C
      LOGICAL ABEL,LOCFRAME,GOPARR,DSKWRK,MASWRK
      CHARACTER*3 AGROUP
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SYMMOL/ GROUP,COMPLEX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     >                T(432),INVT(48),NT
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /FRAME / U1,U2,U3,V1,V2,V3,W1,W2,W3,X0,Y0,Z0
C     COMMON /FRAMES/ X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,DIRECT
C
      DIMENSION GRP(19)
C
      DATA GRP /5HC1   ,5HCS   ,5HCI   ,5HCN   ,5HS2N  ,5HCNH  ,
     >          5HCNV  ,5HDN   ,5HDNH  ,5HDND  ,5HCINFV,5HDINFH,5HT    ,
     >          5HTH   ,5HTD   ,5HO    ,5HOH   ,5HI    ,5HIH   /
C
C     THIS SUBROUTINE SETS SYMMETRY FACTORS FOR THE OCTANTS
C     *****************************************************
C
C     NON-ABELIAN GROUPS ARE HANDLED AS THE HIGHEST ABELIAN SUBGROUPS
C     WHICH HAVE THE SAME MASTER FRAME. E.G. D3 -> C1, NOT C2.
C     PATHETIC AS IT IS, SOME NON-ABELIAN GROUPS ARE DEFILED TO C1.
C
C     CHECK IF THE USER TRIED TO FIDDLE WITH THE LOCAL FRAME.
C     IF THAT IS CORRECT THEN WE CANNOT USE OCTANT SYMMETRY SINCE
C     THEN IT MUST BE REDEFINED APPROPRIATELY (NOT IMPOSSIBLE THOUGH).
C
      IF(ABS(U1-ONE).LT.TM6.AND.ABS(U2).LT.TM6.AND.ABS(U3).LT.TM6.AND.
     *   ABS(V1).LT.TM6.AND.ABS(V2-ONE).LT.TM6.AND.ABS(V3).LT.TM6.AND.
     *   ABS(W1).LT.TM6.AND.ABS(W2).LT.TM6.AND.ABS(W3-ONE).LT.TM6) THEN
         LOCFRAME=.FALSE.
      ELSE
         LOCFRAME=.TRUE.
         IF(MASWRK) WRITE(IW,1010)
      ENDIF
C***********************************************************************
C
C     READ ATOMIC COORDINATES AS 'COORD'
C
C***********************************************************************
      ISYMXY=0
      ISYMXZ=0
      ISYMYZ=0
      ISYMRX=0
      ISYMRY=0
      ISYMRZ=0
      ISYMI =0
C
C     ----- C1 SYMMETRY ; READ ONLY IIII -----
C           CN  FOR ODD N
C           CNV FOR ODD N
C           DN  FOR ODD N
C
      ISYMC1=0
      IF(GROUP.EQ.GRP(1).OR.NT.EQ.1.OR.LOCFRAME.OR.
     *   GROUP.EQ.GRP(4).AND.MOD(NAXIS,2).EQ.1.OR.
     *   GROUP.EQ.GRP(7).AND.MOD(NAXIS,2).EQ.1.OR.
     *   GROUP.EQ.GRP(8).AND.MOD(NAXIS,2).EQ.1) THEN
         AGROUP='C1'
         ISYMC1=1
         LTHE=1
         LPHI=1
C
C     ISYMC1 IS USED ONLY FOR THE OCTANT SYMMETRY.
C
C     ----- D2H SYMMETRY -----
C           DNH FOR EVEN N
C           TH
C           OH
C
      ELSEIF(GROUP.EQ.GRP(9).AND.NAXIS.EQ.2.OR.
     *       GROUP.EQ.GRP(9).AND.MOD(NAXIS,2).EQ.0.OR.
     *       GROUP.EQ.GRP(14).OR.
     *       GROUP.EQ.GRP(17))THEN
         AGROUP='D2H'
         ISYMXY=1
         ISYMXZ=1
         ISYMYZ=1
         ISYMRX=1
         ISYMRY=1
         ISYMRZ=1
         ISYMI =1
         LTHE=2
         LPHI=4
C
C     ----- C2H SYMMETRY -----
C           CNH FOR EVEN N
C
      ELSEIF(GROUP.EQ.GRP(6).AND.NAXIS.EQ.2.OR.
     *       GROUP.EQ.GRP(6).AND.MOD(NAXIS,2).EQ.0) THEN
         AGROUP='C2H'
         ISYMXY=1
         ISYMRZ=1
         ISYMI =1
         LTHE=2
         LPHI=4
C        SINCE PHI IS DIVIDED INTO 4 PARTS AND ONLY TWO ARE USED.
C
C     ----- C2V SYMMETRY ------
C           CNV FOR EVEN N
C
      ELSEIF(GROUP.EQ.GRP(7).AND.NAXIS.EQ.2.OR.
     *       GROUP.EQ.GRP(7).AND.MOD(NAXIS,2).EQ.0) THEN
         AGROUP='C2V'
         ISYMXZ=1
         ISYMYZ=1
         ISYMRZ=1
         LTHE=1
         LPHI=4
C
C     ----- D2 SYMMETRY -----
C           DN FOR EVEN N
C           D2D FOR EVEN N
C           T
C           TD
C           O
C
      ELSEIF(GROUP.EQ.GRP(8).AND.NAXIS.EQ.2.OR.
     *       GROUP.EQ.GRP(8).AND.MOD(NAXIS,2).EQ.0.OR.
     *       GROUP.EQ.GRP(10).AND.MOD(NAXIS,2).EQ.0.OR.
     *       GROUP.EQ.GRP(13).OR.
     *       GROUP.EQ.GRP(15).OR.
     *       GROUP.EQ.GRP(16))THEN
         AGROUP='D2'
         ISYMRX=1
         ISYMRY=1
         ISYMRZ=1
         LTHE=2
         LPHI=2
C
C     ----- C2 SYMMETRY -----
C           CN FOR EVEN N
C           S2N FOR EVEN N
C
      ELSEIF(GROUP.EQ.GRP(4).AND.NAXIS.EQ.2.OR.
     *       GROUP.EQ.GRP(4).AND.MOD(NAXIS,2).EQ.0.OR.
     *       GROUP.EQ.GRP(5).AND.MOD(NAXIS,2).EQ.0) THEN
         AGROUP='C2'
         ISYMRZ=1
         LTHE=1
         LPHI=2
C
C     ----- CI SYMMETRY -----
C           S2N FOR ODD N
C           DND FOR ODD N
C
      ELSEIF(GROUP.EQ.GRP(3).OR.
     *       GROUP.EQ.GRP(5).AND.MOD(NAXIS,2).EQ.1.OR.
     *       GROUP.EQ.GRP(10).AND.MOD(NAXIS,2).EQ.1) THEN
         AGROUP='CI'
         ISYMI =1
         LTHE=1
         LPHI=4
C
C     ----- CS SYMMETRY -----
C           CNH FOR ODD N
C           DNH FOR ODD N
C
      ELSEIF(GROUP.EQ.GRP(2).OR.
     *       GROUP.EQ.GRP(6).AND.MOD(NAXIS,2).EQ.1.OR.
     *       GROUP.EQ.GRP(9).AND.MOD(NAXIS,2).EQ.1) THEN
         AGROUP='CS'
         ISYMXY=1
         LTHE=2
         LPHI=1
      ELSE
        IF(MASWRK) WRITE(IW,1007) GROUP,NAXIS
        CALL ABRT
      ENDIF
C     CHECK IF THE OCTANT SYMMETRY IS CONSISTENT WITH NTHE AND NPHI.
C     THIS IS BECAUSE WE DIVIDE NTHE AND/OR NPHI BY 2 OR 4.
      IF(MOD(NTHE,LTHE).NE.0.OR.MOD(NPHI,LPHI).NE.0.OR.
     *   MOD(NTHE0,LTHE).NE.0.OR.MOD(NPHI0,LPHI).NE.0) THEN
        IF(MASWRK) WRITE(IW,1000) LTHE,LPHI
        ISYMC1=1
      ENDIF
      RETURN
 1000 FORMAT(1X,'FOR GREATER EFFICIENCY, NTHE HAS TO BE A MULTIPLE OF ',
     *       I1,' AND NPHI OF ',I1/1X,'SIMILARLY, NTHE0 AND NPHI0.')
 1007 FORMAT(/'CANNOT USE POINT GROUP ',A8,I2)
 1010 FORMAT(//5X,'WARNING: LOCAL FRAME CHANGE IS DETECTED!',
     *       /5X,'THIS IS LIKELY TO MAKE THE DFT CODE SLOWER.',/)
      END
C*MODULE DFTGRD  *DECK RADPT
      SUBROUTINE RADPT(PTRAD,WTRAD,NRAD)
C**********************************************************************
C
C     DETERMINE THE RADIAL POINTS AND WEIGHTS
C
C**********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER ( ONE=1.0D+00,TWO=2.0D+00)
      DIMENSION PTRAD(NRAD),WTRAD(NRAD)
C
      XTMP=NRAD
      XNP1=XTMP+ONE
      DO 10 IRAD=1,NRAD
         XI=IRAD
C
C     ----- DETERMINE RADIAL POINT -----
C
         PTRAD(IRAD)=(XI*XI)/((XNP1-XI)*(XNP1-XI))
C
C     ----- CALCULATE WEIGHT FOR RADIAL QUADRATURE -----
C
         DEN=(XNP1-XI)*(XNP1-XI)*(XNP1-XI)*(XNP1-XI)
     >        *(XNP1-XI)*(XNP1-XI)*(XNP1-XI)
         WTRAD(IRAD)=(TWO*XNP1*XI*XI*XI*XI*XI)/DEN
C
 10   CONTINUE
      RETURN
      END
C*MODULE DFTGRD  *DECK SYMFCT
      SUBROUTINE SYMFCT(ISYMYZ,ISYMXZ,ISYMXY,ISYMRX,ISYMRY,
     >                  ISYMRZ,ISYMI,NCNTR,
     >                  IQ1,IQ2,IQ3,IQ4,IQ5,IQ6,IQ7,IQ8,
     >                  IFACT,ISYMC1)
C***********************************************************************
C
C     THIS ROUTINE DETERMINES WHICH QUADRANTS TO INTEGRATE OVER AND THE
C     FACTOR TO MULTIPLY THE SYMMETRY INTEGRATION BY.
C
C     ATOMS THAT LIE IN SYMMETRY PLANES OR ARE ON THE SYMMETRY AXES
C     (IN THE ABELIAN SUBGROUP) CAN HAVE ONLY A FRACTION OF THE GRID PTS
C     EXPLICITLY COMPUTED AND THE REST DUPLICATED WITH SYMMETRY FACTORS.
C     THIS IS CALLED "OCTANT SYMMETRY" AND IS TREATED HERE.
C     ATOMS THAT ARE CONNECTED TO A SYMMETRY UNIQUE ATOM BY SYMMETRY
C     OPERATIONS (IN THE FULL GROUP) ARE ALSO NOT COMPUTED AND THE GRID
C     INTEGRATION OVER A SYMMETRY UNIQUE ATOM IS MULTIPLIED BY A FACTOR.
C     THIS IS HANDLED IN SYMUNQ.
C
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, THRSH=1.0D-10)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),COORD(3,MXATM),IAN(MXATM)
C
      DIMENSION IFACT(NAT)
C***********************************************************************
C
C     QUADRANTS ARE LABELED
C     XYZ=1, XY(-Z)=2, X(-Y)Z=3, X(-Y)(-Z)=4
C     (-X)YX=5, (-X)Y(-Z)=6, (-X)(-Y)Z=7, (-X)(-Y)(-Z)=8
C
C     DETERMINE WHICH SYMMETRY ELEMENTS APPLY TO SYMMETRY-UNIQUE
C     ATOM NCNTR
C
C***********************************************************************
      IYZ=ISYMYZ
      IXZ=ISYMXZ
      IXY=ISYMXY
      IRX=ISYMRX
      IRY=ISYMRY
      IRZ=ISYMRZ
      IMI=ISYMI
C***********************************************************************
C     C1 SYMMETRY
C***********************************************************************
      IF(ISYMC1.EQ.1) THEN
        IFACT(NCNTR) = 1
        GOTO 300
      ENDIF
C***********************************************************************
C
C     SET SYMMETRICAL FACTORS
C
C***********************************************************************
      IF(ABS(COORD(1,NCNTR)).GT.THRSH.AND.
     >   ABS(COORD(2,NCNTR)).GT.THRSH.AND.
     >   ABS(COORD(3,NCNTR)).GT.THRSH)THEN
         ISYMYZ=0
         ISYMXZ=0
         ISYMXY=0
         ISYMRX=0
         ISYMRY=0
         ISYMRZ=0
         ISYMI=0
C***********************************************************************
      ELSE IF(ABS(COORD(1,NCNTR)).LE.THRSH.AND.
     >   ABS(COORD(2,NCNTR)).GT.THRSH.AND.
     >   ABS(COORD(3,NCNTR)).GT.THRSH)THEN
         ISYMXZ=0
         ISYMXY=0
         ISYMRX=0
         ISYMRY=0
         ISYMRZ=0
         ISYMI=0
C***********************************************************************
      ELSE IF(ABS(COORD(1,NCNTR)).GT.THRSH.AND.
     >   ABS(COORD(2,NCNTR)).LE.THRSH.AND.
     >   ABS(COORD(3,NCNTR)).GT.THRSH)THEN
         ISYMYZ=0
         ISYMXY=0
         ISYMRX=0
         ISYMRY=0
         ISYMRZ=0
         ISYMI=0
C***********************************************************************
      ELSE IF(ABS(COORD(1,NCNTR)).GT.THRSH.AND.
     >   ABS(COORD(2,NCNTR)).GT.THRSH.AND.
     >   ABS(COORD(3,NCNTR)).LE.THRSH)THEN
         ISYMYZ=0
         ISYMXZ=0
         ISYMRX=0
         ISYMRY=0
         ISYMRZ=0
         ISYMI=0
C***********************************************************************
      ELSE IF(ABS(COORD(1,NCNTR)).LE.THRSH.AND.
     >   ABS(COORD(2,NCNTR)).LE.THRSH.AND.
     >   ABS(COORD(3,NCNTR)).GT.THRSH)THEN
         ISYMXY=0
         ISYMRX=0
         ISYMRY=0
         ISYMI=0
C***********************************************************************
      ELSE IF(ABS(COORD(1,NCNTR)).LE.THRSH.AND.
     >   ABS(COORD(2,NCNTR)).GT.THRSH.AND.
     >   ABS(COORD(3,NCNTR)).LE.THRSH)THEN
         ISYMXZ=0
         ISYMRX=0
         ISYMRZ=0
         ISYMI=0
C***********************************************************************
      ELSE IF(ABS(COORD(1,NCNTR)).GT.THRSH.AND.
     >   ABS(COORD(2,NCNTR)).LE.THRSH.AND.
     >   ABS(COORD(3,NCNTR)).LE.THRSH)THEN
         ISYMYZ=0
         ISYMRY=0
         ISYMRZ=0
         ISYMI=0
      ENDIF
C***********************************************************************
C
C     INITIAL SYMMETRY FACTOR
C
C***********************************************************************
      IFACT(NCNTR)=1
C***********************************************************************
C
C     DETERMINE WHICH QUADRANTS ARE EQUIVALENT
C
C***********************************************************************
      IQ1=1
      IQ2=1
      IQ3=1
      IQ4=1
      IQ5=1
      IQ6=1
      IQ7=1
      IQ8=1
C***********************************************************************
      IF(ISYMYZ.EQ.1.OR.ISYMRY.EQ.1.OR.ISYMRZ.EQ.1.OR.ISYMI.EQ.1)THEN
         IQ5=0
         IQ6=0
         IQ7=0
         IQ8=0
         IF(ISYMXZ.EQ.1.OR.ISYMRX.EQ.1)THEN
            IQ3=0
            IQ4=0
            IF(ISYMXY.EQ.1)THEN
               IQ2=0
               IFACT(NCNTR)=8
               GOTO 300
            ENDIF
         IFACT(NCNTR)=4
         GO TO 300
         ENDIF
         IF(ISYMXY.EQ.1)THEN
            IQ2=0
            IQ4=0
            IFACT(NCNTR)=4
            GOTO 300
         ENDIF
         IFACT(NCNTR)=2
      ENDIF
      IF(ISYMXZ.EQ.1.OR.ISYMRX.EQ.1)THEN
            IQ3=0
            IQ4=0
            IQ7=0
            IQ8=0
            IF(ISYMXY.EQ.1)THEN
               IQ2=0
               IQ6=0
               IFACT(NCNTR)=4
               GOTO 300
            ENDIF
      IFACT(NCNTR)=2
      GOTO 300
      ENDIF
      IF(ISYMXY.EQ.1)THEN
         IQ2=0
         IQ4=0
         IQ6=0
         IQ8=0
         IFACT(NCNTR)=2
      ENDIF
  300 CONTINUE
      ISYMYZ=IYZ
      ISYMXZ=IXZ
      ISYMXY=IXY
      ISYMRX=IRX
      ISYMRY=IRY
      ISYMRZ=IRZ
      ISYMI=IMI
      RETURN
      END
C*MODULE DFTGRD  *DECK SYMUNQ
      SUBROUTINE SYMUNQ(IFACTR,INATM,NSYMAT,NDEG,NEQATM,IUNIQUE,IZCORE)
C***********************************************************************
C
C      THIS ROUTINE SETS UP THE FACTOR
C         FOR THE SYMMETRY-UNIQUE ATOM TO BE MULTIPLIED BY.
C
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000,MXATM=2000)
C
      DIMENSION IFACTR(NAT),INATM(NAT),NDEG(NAT),NEQATM(NAT,NAT),
     *          IUNIQUE(NAT),IZCORE(*)
C
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
C     DETERMINE SYMMETRY UNIQUE ATOMS AND SYMMETRY FACTORS
C
      IF (NT.EQ.1) THEN
         DO 5 I=1,NAT
            IFACTR(I)=1
            INATM(I)=I
            IUNIQUE(I)=I
            INUC=INT(ZAN(I)+IZCORE(I))
            IF(NFG.NE.0) INUC=IAN(I)
            IF(INUC.EQ.0) IUNIQUE(I)=0
 5       CONTINUE
         NSYMAT=NAT
         RETURN
      ENDIF
C
      NSYMAT=0
      CALL VICLR(IFACTR,1,NAT)
      CALL VICLR(INATM,1,NAT)
      CALL VICLR(NDEG,1,NAT)
      CALL VICLR(NEQATM,1,NAT*NAT)
      CALL VICLR(IUNIQUE,1,NAT)
      JATM=1
 10   CONTINUE
      IF (JATM.GT.NAT) THEN
         GOTO 90
      ELSEIF (JATM.NE.1) THEN
         DO 20 KATM=1,JATM-1
            IF (NDEG(KATM).EQ.0) GOTO 20
            DO 30 IDEG=1,NDEG(KATM)
               IF (JATM.EQ.NEQATM(KATM,IDEG)) THEN
                  JATM=JATM+1
                  GOTO 10
               ENDIF
 30         CONTINUE
 20      CONTINUE
      ENDIF
C
      NSYMAT=NSYMAT+1
      NDEG(JATM)=1
      INATM(NSYMAT)=JATM
C
C
      DO 40 IOPER=1,NT-1
         DO 50 JOPER=1,IOPER
           IF (MAPCTR(JATM,JOPER).EQ.MAPCTR(JATM,IOPER+1)) GOTO 40
 50      CONTINUE
         NEQATM(JATM,NDEG(JATM))=MAPCTR(JATM,IOPER+1)
         NDEG(JATM)=NDEG(JATM)+1
C
 40   CONTINUE
      IFACTR(NSYMAT)=NDEG(JATM)
C
      JATM=JATM+1
      GOTO 10
C
 90   CONTINUE
      DO INC = 1, NSYMAT
        IUNIQUE(INATM(INC))=INC
        I=INATM(INC)
        INUC=INT(ZAN(I)+IZCORE(I))
        IF(NFG.NE.0) INUC=IAN(I)
        IF(INUC.EQ.0) IUNIQUE(I)=0
      ENDDO
      RETURN
      END
C*MODULE DFTGRD  *DECK DFTTRF
      SUBROUTINE DFTTRF(UROHF,L1,NA,NB,CA,CB,AOX,MOXA,MOXB,TOL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION MOXA,MOXB
      LOGICAL UROHF
      DIMENSION CA(L1,NA),CB(L1,NB),AOX(L1),MOXA(L1),MOXB(L1)
C
C     AO->MO TRANSFORMATION
C     (FOR SPEED, RUN ONE RATHER THAN TWO LOOPS, SO ZERO EXTRA -MOXB-)
C     ALL THIS CAN BE COMPRESSED INTO A SINGLE LINE:
C     CALL MRARBR(AOX,1,1,L1,COEFFB,L1,NOB,VMOB,1)
C
                CALL VCLR(MOXA,1,NA)
      IF(UROHF) CALL VCLR(MOXB,1,NA)
C
      DO I=1,L1
        AO=AOX(I)
        IF(ABS(AO).GT.TOL) THEN
          IF(UROHF) THEN
             DO J=1,NA
                MOXA(J)=MOXA(J)+AO*CA(I,J)
                MOXB(J)=MOXB(J)+AO*CB(I,J)
             ENDDO
          ELSE
             DO J=1,NA
                MOXA(J)=MOXA(J)+AO*CA(I,J)
             ENDDO
          ENDIF
        ENDIF
      ENDDO
C--   IF(.NOT.UROHF) CALL DCOPY(NB,MOXA,1,MOXB,1)
      RETURN
      END
C*MODULE DFTGRD  *DECK DFTTRFG
      SUBROUTINE DFTTRFG(UROHF,L1,NA,NB,CA,CB,
     *                   AOX,GAOX,GAOY,GAOZ,GMOXA,
     *                   GMOYA,GMOZA,GMOXB,GMOYB,GMOZB,TOL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL UROHF
      DIMENSION CA(L1,NA),CB(L1,NB),AOX(L1),GAOX(L1),GAOY(L1),GAOZ(L1),
     *      GMOXA(L1),GMOYA(L1),GMOZA(L1),GMOXB(L1),GMOYB(L1),GMOZB(L1)
C
C     AO->MO TRANSFORMATION FOR THE GRADIENT TERMS.
C     SIMILARLY, A SINGLE CALL TO MRARBR WOULD DO THE TRICK.
C     THIS ROUTINE ASSUMES THAT CB (BETA-LCAO COEFFICIENTS) ARE ZERO
C     FOR "MISSING" ORBITALS, THAT IS, LAST NA-NB ORBITALS, FOR UROHF.
C
      CALL VCLR(GMOXA,1,NA)
      CALL VCLR(GMOYA,1,NA)
      CALL VCLR(GMOZA,1,NA)
C        SHOULD BE NB IN PRINCIPLE
      IF(UROHF) THEN
         CALL VCLR(GMOXB,1,NA)
         CALL VCLR(GMOYB,1,NA)
         CALL VCLR(GMOZB,1,NA)
      ENDIF
C
C--   TOL2=TOL*TOL
C
      DO I=1,L1
        IF(ABS(AOX(I)).GT.TOL) THEN
C         NN=NC(I)
        GX=GAOX(I)
        GY=GAOY(I)
        GZ=GAOZ(I)
C       IF(ABS(GX)+ABS(GY)+ABS(GZ).GT.TOL) THEN
          IF(UROHF) THEN
          DO J=1,NA
            COA=CA(I,J)
            COB=CB(I,J)
            GMOXA(J)=GMOXA(J)+GX*COA
            GMOYA(J)=GMOYA(J)+GY*COA
            GMOZA(J)=GMOZA(J)+GZ*COA
            GMOXB(J)=GMOXB(J)+GX*COB
            GMOYB(J)=GMOYB(J)+GY*COB
            GMOZB(J)=GMOZB(J)+GZ*COB
          ENDDO
          ELSE
          DO J=1,NA
            COA=CA(I,J)
            GMOXA(J)=GMOXA(J)+GX*COA
            GMOYA(J)=GMOYA(J)+GY*COA
            GMOZA(J)=GMOZA(J)+GZ*COA
          ENDDO
          ENDIF
        ENDIF
      ENDDO
C
CNB   THIS PART CAN PROBABLY BE NEGLECTED
C     IF(.NOT.UROHF) THEN
C       CALL DCOPY(NB,GMOXA,1,GMOXB,1)
C       CALL DCOPY(NB,GMOXA,1,GMOXB,1)
C       CALL DCOPY(NB,GMOXA,1,GMOXB,1)
C     ENDIF
      RETURN
      END
C*MODULE DFTGRD  *DECK SWGRID
      SUBROUTINE SWGRID
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0
C
C     SWITCH THE FINE AND COARSE GRID SIZES
C
      NRADT=NRAD
      NTHET=NTHE
      NPHIT=NPHI
      NRAD=NRAD0
      NTHE=NTHE0
      NPHI=NPHI0
      NRAD0=NRADT
      NTHE0=NTHET
      NPHI0=NPHIT
C
      RETURN
      END
C*MODULE DFTGRD  *DECK DEREXC
      SUBROUTINE DEREXC(TOTWT,NAPTS,COEFFA,COEFFB,DA,DB,IANGN,IFACTR,
     *                  INC0,VMOA,DVMOXA,DVMOYA,DVMOZA,VMOB,DVMOXB,
     *                  DVMOYB,DVMOZB,WGHT,RI,RIJ,WTINTR,AIJ,WTRAD,
     *                  NCNTR,L1,TOTELE,NANG,PTRAD,XDAT,YDAT,ZDAT,
     *                  ATMXVEC,ATMYVEC,ATMZVEC,ANGXVL,ANGYVL,ANGZVL,
     *                  BSLRDII,RSQRD,PCOEFF,EXPS,AOX,GAOX,GAOY,GAOZ,
     *                  G2AOXX,G2AOYY,G2AOZZ,G2AOXY,G2AOYZ,G2AOXZ,
     *                  IJX,IJY,IJZ,BOHR,NEEDGR,DEDFT,IZCORE,UVEC,
     *                  DWTINT,DWTTOT,DTOTWT,OUT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL UROHF,OUT,NEEDGR,GOPARR,DSKWRK,MASWRK,NXT
C
      PARAMETER (MXATM=2000,ZERO=0.0D+00,ONE=1.0D+00,TWO=2.0D+00)
C
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DIMENSION TOTWT(NAT,*),COEFFA(L1,*),COEFFB(L1,*),
     >          VMOA(NUM),DVMOXA(NUM),DVMOYA(NUM),DVMOZA(NUM),
     >          VMOB(NUM),DVMOXB(NUM),DVMOYB(NUM),DVMOZB(NUM),
     *          DA(NUM,NUM),DB(NUM,NUM)
      DIMENSION NAPTS(NAT),IANGN(NAT,2),IFACTR(NAT)
      DIMENSION AOX(NUM),GAOX(NUM),GAOY(NUM),GAOZ(NUM),
     *          G2AOXX(NUM),G2AOYY(NUM),G2AOZZ(NUM),
     *          G2AOXY(NUM),G2AOYZ(NUM),G2AOXZ(NUM),
     *          DEDFT(3,NAT),EXPS(*),IZCORE(*),
     *          PTRAD(*),XDAT(NAT,NTHE*NPHI),YDAT(NAT,NTHE*NPHI),
     *          ZDAT(NAT,NTHE*NPHI),RSQRD(*),PCOEFF(*),BSLRDII(*),
     *          ATMXVEC(NAT,NAT),ATMYVEC(NAT,NAT),ATMZVEC(NAT,NAT),
     *          ANGXVL(NAT,*),ANGYVL(NAT,*),ANGZVL(NAT,*),
     *          WGHT(NAT,NTHE*NPHI),RI(NAT),RIJ(NAT,NAT),WTINTR(NAT),
     *          AIJ(NAT,NAT),WTRAD(*),UVEC(3,NAT),DWTINT(3,NAT,NAT),
     *          DWTTOT(3,NAT),DTOTWT(3,NAT)
      DIMENSION IJX(35),IJY(35),IJZ(35)
C
      DATA DEBUG/8HDEBUG   /, DFTDER/8HDFTDER  /
      DATA UHF,ROHF/8HUHF     ,8HROHF    /
C
C     THIS ROUTINE IS A CLONE OF DMATD, SEE ADDITIONAL REFERENCES THERE.
C     HERE A CONTRIBUTION TO THE ENERGY GRADIENT DUE TO INTEGRATION OVER
C     GRID POINTS SURROUNDING ATOM NCNTR IS COMPUTED.
C
      UROHF = SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
      OUT = EXETYP.EQ.DFTDER  .OR.  EXETYP.EQ.DEBUG
C
      TOTEXC=ZERO
C
      NUC=INT(ZAN(NCNTR)+IZCORE(NCNTR))
      IF(NFG.NE.0) NUC=IAN(NCNTR)
      RAD=BSLRDII(NUC)/BOHR
      EXEC = ZERO
      NOA = NA
      NOB = NB
      NPT=NRAD*NTHE*NPHI
C
C     SET CUT-OFFS FOR DFT GRADIENT:
C     WCUTOFF - CELL VOLUME (WEIGHT)
C     RCUTOFF - ELECTRON DENSITY
C     CCUTOFF - LCAO COEFFICIENTS (AS USED FOR AO -> MO TRANSFORMATION)
C     FCUTOFF - FUNCTIONAL DERIVATIVE
C     AT PRESENT WEIGHT DERIVATIVE CUTOFF IS FOUND UNUSEFUL.
C
C     DFTTHR IS NOT USED IN THE GRADIENT.
C     INSTEAD, DFTGTHR IS USED AS FOLLOWS:
C     DFTGTHR=1 USE DEFAULT THRESHOLDS
C     DFTGTHR<1 SET ALL THRESHOLDS TO DFTGTHR
C     DFTGTHR>1 DIVIDE DEFAULT THRESHOLDS BY DFTGTHR (MAKE TIGHTER).
C     DFTTHR SET TO 1E-15 AS BEFORE SETS ALL GRAD THRESHOLDS TO 1E-15.
C
C     IF HIGHER ACCURACY IN GRADIENT DUE TO THRESHOLDS IS THOUGHT THE
C     RECOMMENDED WAY IS TO SET DFTGTHR=10 (WHICH MEANS ADD ABOUT
C     1 DIGIT OF ACCURACY, RELATIVE TO NO THRESHOLDS; DFT GRID
C     INTEGRATION HAS ITS OWN ERRORS, FINALLY, 2E INTEGRALS ALSO HAVE
C     CERTAIN ERROR BARS ON THEM).
C
      THRV1=ZERO
      THRV2=1.0D+00
      IF(DFTGTHR.GT.1.01D+00) THEN
        THRV2=DFTGTHR
      ELSE IF(DFTGTHR.LT.0.99D+00) THEN
        THRV1=DFTGTHR
      ENDIF
      IF(DFTTHR.LT.1.1D-15.AND.DFTTHR.NE.ZERO) THRV1=1.0D-15
C     THRV2=1.0D+00
      WCUTOFF=3.0D-05/NPT/THRV2
      RCUTOFF=1.0D-02/NPT/THRV2
      CCUTOFF=1.0D-01/NPT/THRV2
      FCUTOFF=3.0D-04/NPT/THRV2
C     WDCUTOFF=1.0D+01*WCUTOFF
C     IF(DFTTHR.LT.1.1D-15.AND.DFTTHR.NE.ZERO) THEN
      IF(THRV1.NE.ZERO) THEN
C        THRV1=1.0D-15
         WCUTOFF=THRV1
         RCUTOFF=THRV1
         CCUTOFF=THRV1
         FCUTOFF=THRV1
C        WCUTOFF=WCUTOFF*1D-1
C        RCUTOFF=RCUTOFF*1D-1
C        CCUTOFF=CCUTOFF*1D-1
C        FCUTOFF=FCUTOFF*1D-1
C        WDCUTOFF=WDCUTOFF*1D-1
      ENDIF
C     WRITE(6,*) 'WWWG',WCUTOFF,RCUTOFF,CCUTOFF,FCUTOFF,DFTGTHR
C
C     ANGXVL=X**I, ANGYVL=Y**J, ANGZVL=Z**K, FOR ALL VALUES OF I,J,K
C     WHERE (X,Y,Z) IS THE CENTRE OF A DFT GRID POINT RI
C     MINUS AN ATOMIC CENTRE RA: XYZ= (RI - RA)
C
      CALL VCLR(ANGXVL(1,1),1,NAT)
      CALL VCLR(ANGYVL(1,1),1,NAT)
      CALL VCLR(ANGZVL(1,1),1,NAT)
      CALL DACOPY(NAT,ONE,ANGXVL(1,2),1)
      CALL DACOPY(NAT,ONE,ANGYVL(1,2),1)
      CALL DACOPY(NAT,ONE,ANGZVL(1,2),1)
C
C     ---- SYMMETRY MULTIPLICATION FACTOR ----
C
C      FACT = IIFACT(NCNTR)*IFACTR(INC0)
      FACT=IFACTR(INC0)
C      FACT=1
C
C     DYNAMIC LOAD BALANCING: DIVIDE ALL POINTS INTO NGRAN*NPROC CHUNKS.
C     NGRAN CANNOT BE SET IN THE INPUT FILE. THEN EACH NODE GETS A CHUNK
C     AND TOILS IT OVER UNTIL READY TO ASK FOR MORE.
C
      LOOP=0
      NGRAN=20
      NLOOP=NRAD*(IANGN(NCNTR,2)-IANGN(NCNTR,1)+1)
      MCHUNK=(NLOOP-1)/(NGRAN*NPROC)+1
      NEXT  = -1
      NXT = IBTYP.EQ.1
C     WRITE(6,*) 'WWWAA',NXT
C***********************************************************************
C     LOOP OVER RADIAL GRIDS
C     LOOP OVER ANGULAR GRIDS
C***********************************************************************
      DO 20 IRADPT = 1, NRAD
C
        RADWT=RAD*RAD*RAD*WTRAD(IRADPT)
        R1= RAD*PTRAD(IRADPT)
        DO 10 IANGPT = IANGN(NCNTR,1), IANGN(NCNTR,2)
          IPT=(IRADPT-1)*NAPTS(NCNTR)+IANGPT
C
          IF(GOPARR) THEN
            IF(NXT) THEN
              LOOP=LOOP+1
              ICHUNK=(LOOP-1)/MCHUNK
              IF(ICHUNK.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
              IF(NEXT.NE.ICHUNK) GOTO 10
            ELSE
              IF(MOD(IPT,NPROC).NE.ME) GOTO 10
            ENDIF
          ENDIF
          IPTME=(IPT-1)/NPROC+1
          IF(NXT) IPTME=IPT
C
          FTOTWT = TOTWT(NCNTR,IPTME)*FACT
          FTOTWT0=FTOTWT/FACT
          IF(FTOTWT0.LT.WCUTOFF) THEN
             GOTO 10
          ENDIF
C
          XD=R1*XDAT(NCNTR,IANGPT)
          YD=R1*YDAT(NCNTR,IANGPT)
          ZD=R1*ZDAT(NCNTR,IANGPT)
          DO 610 IATM=1,NAT
            XCDNT=ATMXVEC(NCNTR,IATM)+XD
            YCDNT=ATMYVEC(NCNTR,IATM)+YD
            ZCDNT=ATMZVEC(NCNTR,IATM)+ZD
            RSQRD(IATM)=XCDNT**2+YCDNT**2+ZCDNT**2
C
C           NORMALISE THE R(A) VECTOR, = R(I)-R(A)
C
            RI(IATM)=SQRT(RSQRD(IATM))
            UVEC(1,IATM) = XCDNT / RI(IATM)
            UVEC(2,IATM) = YCDNT / RI(IATM)
            UVEC(3,IATM) = ZCDNT / RI(IATM)
C
            ANGXVL(IATM,3)=XCDNT
            ANGYVL(IATM,3)=YCDNT
            ANGZVL(IATM,3)=ZCDNT
            DO 35 IANG=3,NANG
               ANGXVL(IATM,IANG+1)=ANGXVL(IATM,IANG)*XCDNT
               ANGYVL(IATM,IANG+1)=ANGYVL(IATM,IANG)*YCDNT
               ANGZVL(IATM,IANG+1)=ANGZVL(IATM,IANG)*ZCDNT
 35         CONTINUE
  610     CONTINUE
C
          CALL DFTAO(IJX,IJY,IJZ,AOX,ANGXVL,ANGYVL,ANGZVL,PCOEFF,EXPS,
     *               RSQRD,NAT,L1)
C***********************************************************************
C     FORM DENSITY AT THIS POINT
C***********************************************************************
          CALL DFTTRF(UROHF,L1,NOA,NOB,COEFFA,COEFFB,AOX,VMOA,VMOB,
     *                CCUTOFF)
          ROA=DDOT(NOA,VMOA,1,VMOA,1)
          IF(UROHF) THEN
            ROB=DDOT(NOB,VMOB,1,VMOB,1)
          ELSE
            ROB=ROA
          ENDIF
          IF(ROA+ROB.LT.RCUTOFF) THEN
             GOTO 10
          ENDIF
C         IF(FTOTWT0.GE.WDCUTOFF) THEN
C
C         COMPUTE THE WEIGHT DERIVATIVE.
C
          WGHTNOW=WGHT(NCNTR,IANGPT)
          NITR=4
          CALL GRDOCT(NAT,NITR,NCNTR,ATMXVEC,ATMYVEC,ATMZVEC,RI,RIJ,AIJ,
     *          WGHTNOW,RADWT,WTINTR,UVEC,DWTINT,DWTTOT,DTOTWT)
C         ELSE
C           CALL VCLR(DTOTWT,1,3*NAT)
C         ENDIF
C
          CALL DFTGAO(IJX,IJY,IJZ,GAOX,GAOY,GAOZ,ANGXVL,ANGYVL,ANGZVL,
     *                PCOEFF,EXPS,RSQRD,NAT,L1)
C         ---------------
          IF(NEEDGR) THEN
C         ----------------
          CALL DFTG2AO(IJX,IJY,IJZ,G2AOXX,G2AOYY,G2AOZZ,G2AOXY,G2AOYZ,
     *          G2AOXZ,ANGXVL,ANGYVL,ANGZVL,PCOEFF,EXPS,RSQRD,NAT,L1)
C***********************************************************************
C     FORM DENSITY GRADIENT AT THIS POINT
C***********************************************************************
            CALL DFTTRFG(UROHF,L1,NOA,NOB,COEFFA,COEFFB,
     *                   AOX,GAOX,GAOY,GAOZ,
     *                   DVMOXA,DVMOYA,DVMOZA,DVMOXB,DVMOYB,DVMOZB,
     *                   CCUTOFF)
            GRADXA=TWO*DDOT(NOA,VMOA,1,DVMOXA,1)
            GRADYA=TWO*DDOT(NOA,VMOA,1,DVMOYA,1)
            GRADZA=TWO*DDOT(NOA,VMOA,1,DVMOZA,1)
            IF(UROHF) THEN
              GRADXB=TWO*DDOT(NOB,VMOB,1,DVMOXB,1)
              GRADYB=TWO*DDOT(NOB,VMOB,1,DVMOYB,1)
              GRADZB=TWO*DDOT(NOB,VMOB,1,DVMOZB,1)
            ELSE
              GRADXB=GRADXA
              GRADYB=GRADYA
              GRADZB=GRADZA
            ENDIF
C***********************************************************************
C      FORM GRADIENT INVARIENT (GRAD DOT GRAD)
C***********************************************************************
            GRDAA  = GRADXA*GRADXA+GRADYA*GRADYA+GRADZA*GRADZA
            GRDBB  = GRADXB*GRADXB+GRADYB*GRADYB+GRADZB*GRADZB
            GRDAB  = GRADXA*GRADXB+GRADYA*GRADYB+GRADZA*GRADZB
          ENDIF
C***********************************************************************
C     THE EXCHANGE-CORRELATION FUNCTIONAL CALCULATION
C***********************************************************************
          VXCA1=ZERO
          VXCB1=ZERO
          DUMAX=ZERO
          DUMAY=ZERO
          DUMAZ=ZERO
          DUMBX=ZERO
          DUMBY=ZERO
          DUMBZ=ZERO
          XALPHA=ZERO
          XGRD=ZERO
          ECF=ZERO
          CALL CALCEXC(ROA,ROB,FTOTWT,GRDAA,GRDBB,GRDAB,
     >                 GRADXA,GRADYA,GRADZA,GRADXB,GRADYB,GRADZB,
     >                 XALPHA,XGRD,VXCA1,DUMAX,DUMAY,DUMAZ,
     >                 VXCB1,DUMBX,DUMBY,DUMBZ,ECF)
          DUMA = VXCA1
          DUMB = VXCB1
          EXEC1= XALPHA + XGRD + ECF
          EXEC = EXEC + EXEC1
          IF(NOB.EQ.0) THEN
             DUMB=ZERO
             GRADXB=ZERO
             GRADYB=ZERO
             GRADZB=ZERO
          ENDIF
          UPBOW=ABS(DUMA)*TWO
          IF(UROHF) UPBOW=ABS(DUMA)+ABS(DUMB)
C         IF(NEEDGR) THEN
C           UPBOW=UPBOW+ABS(DUMAX)+ABS(DUMAY)+ABS(DUMAZ)
C           IF(UROHF) UPBOW=UPBOW+ABS(DUMBX)+ABS(DUMBY)+ABS(DUMBZ)
C         ENDIF
          UPBOW=UPBOW*TWO*FTOTWT0
          IF(UPBOW.GT.FCUTOFF) THEN
C         XGRD1=XGRD1+XGRD
C***********************************************************************
C      COMPUTE DFT GRADIENT
C***********************************************************************
          CALL DFTGDV(NCNTR,DEDFT,FTOTWT,NUM,DUMA,DUMAX,DUMAY,DUMAZ,
     >         DUMB,DUMBX,DUMBY,DUMBZ,AOX,GAOX,GAOY,GAOZ,
     >         G2AOXX,G2AOYY,G2AOZZ,G2AOXY,G2AOYZ,G2AOXZ,
     >         DA,DB,NEEDGR,UROHF)
          ENDIF
C
C         ADD THE WEIGHT DERIVATIVE CONTRIBUTION
C         (WEIGHT DERIVATIVE ITSELF TIMES THE FUNCTIONAL VALUE).
C
          EXEC2=EXEC1/FTOTWT*FACT
C         WRITE(6,*) 'WWWWEIGHT',EXEC2,DTOTWT(1,1),DEDFT(1,1)
          CALL DAXPY(3*NAT,EXEC2,DTOTWT,1,DEDFT,1)
C
C     ----- THE TOTAL ELECTRON DENSITY -----
C
          TOTELE=TOTELE+FTOTWT*(ROA+ROB)
C
   10   CONTINUE
C
C     ----- NEXT RADIAL POINT -----
C
   20 CONTINUE
C
      TOTEXC = TOTEXC+EXEC
C
      IF(GOPARR.AND.NXT) CALL DDI_DLBRESET
      IF (OUT) THEN
         WRITE(IW,9999)NCNTR,TOTEXC,TOTELE
         DO INAT=1,NAT
            WRITE(IW,9998) DEDFT(1,INAT),DEDFT(2,INAT),DEDFT(3,INAT)
         ENDDO
      ENDIF
      RETURN
 9998 FORMAT(/F20.10,F20.10,F20.10)
 9999 FORMAT(/5X,'ATM',I8,'EXC=',F20.10,5X,'TOTELE=',F20.10)
      END
C*MODULE DFTGRD  *DECK OCT
      SUBROUTINE OCT(ITMP,XDAT,YDAT,ZDAT,ATMXVEC,ATMYVEC,ATMZVEC,RI,
     *             RIJ,NAT,AIJ,WTINTR,TOTWT,WGHT,IANGN,NITR,NCNTR,RAD,
     *               NAPTS,PTRAD,WTRAD)
C**********************************************************************
C
C     THIS ROUTINE SETS UP THE NUMERICAL GRID FOR THE XYZ OCTANT
C
C**********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER ( ZERO=0.0D+00, ONE=1.0D+00)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION WGHT(NAT,NTHE*NPHI),XDAT(NAT,NTHE*NPHI),
     >          YDAT(NAT,NTHE*NPHI),ZDAT(NAT,NTHE*NPHI),
     >          ATMXVEC(NAT,NAT),ATMYVEC(NAT,NAT),ATMZVEC(NAT,NAT),
     *          RI(NAT),RIJ(NAT,NAT),WTINTR(NAT),TOTWT(NAT,*),
     *          AIJ(NAT,NAT),IANGN(NAT,2),NAPTS(NAT),PTRAD(ITMP),
     *          WTRAD(ITMP)
C
C**********************************************************************
C
C     DETERMINE GRID POINT VOLUMES (=WEIGHTS) FOR ONE OCTANT.
C     XCDNT ETC CONTAIN COORDINATES OF THE CENTRE OF A GRID POINT
C
C**********************************************************************
C
      NXT = IBTYP.EQ.1
      DO 600 IRAD=1,ITMP
         RADWT=RAD*RAD*RAD*WTRAD(IRAD)
         R1= RAD*PTRAD(IRAD)
C
         DO 605 I=IANGN(NCNTR,1),IANGN(NCNTR,2)
            IPT=(IRAD-1)*NAPTS(NCNTR)+I
C
            IF(MOD(IPT,NPROC).NE.ME) GOTO 605
            IPTME=(IPT-1)/NPROC+1
            IF(NXT) IPTME=IPT
C
            XD=R1*XDAT(NCNTR,I)
            YD=R1*YDAT(NCNTR,I)
            ZD=R1*ZDAT(NCNTR,I)
         DO 610 IATM=1,NAT
            XCDNT=ATMXVEC(NCNTR,IATM)+XD
            YCDNT=ATMYVEC(NCNTR,IATM)+YD
            ZCDNT=ATMZVEC(NCNTR,IATM)+ZD
            RSQRD=XCDNT**2+YCDNT**2+ZCDNT**2
            RI(IATM)=SQRT(RSQRD)
  610    CONTINUE
         WTTOT=ZERO
C
         CALL DACOPY(NAT,ONE,WTINTR,1)
C
         DO 611 IATM=1,NAT
C            WTINTR(IATM)=ONE
            IF (ABS(AIJ(IATM,1)+1.0D+00).LT.1.0D-05) THEN
              WTINTR(IATM) = 0.0D+00
              GOTO 614
            ENDIF
            DO 612 JATM=1,IATM-1
              IF (ABS(AIJ(IATM,JATM)-1.0D+00).LT.1.0D-05) THEN
                CUTIJ = 1.0D+00
                GOTO 620
              ENDIF
              ZMUIJ=(RI(IATM)-RI(JATM))/RIJ(IATM,JATM)
              XMUIJ=ZMUIJ
     >           +AIJ(IATM,JATM)*(ONE-ZMUIJ*ZMUIJ)
C***********************************************************************
C
C     BECKE'S FUZZY CELL METHOD FOR MOLECULAR GRID QUADRATURE
C
C***********************************************************************
C             F4=FUZZY(NITR,XMUIJ)
              F4=XMUIJ
              DO K=1,NITR
                F4=F4*(1.5D+00-0.5D+00*F4*F4)
              END DO
              F2=0.5D+00*F4
              CUTIJ=0.5D+00-F2
              CUTJI=0.5D+00+F2
  620       CONTINUE
              WTINTR(IATM)=WTINTR(IATM)*CUTIJ
              WTINTR(JATM)=WTINTR(JATM)*CUTJI
  612       CONTINUE
  614       CONTINUE
  611    CONTINUE
         DO IATM=1,NAT
            WTTOT=WTTOT+WTINTR(IATM)
         ENDDO
C
         ATMWT=WTINTR(NCNTR)/WTTOT
         TOTWT(NCNTR,IPTME)=ATMWT*RADWT*WGHT(NCNTR,I)
C
  605    CONTINUE
  600 CONTINUE
      RETURN
      END
C*MODULE DFTGRD  *DECK DFTGDV
      SUBROUTINE DFTGDV(NCNTR,DEDFT,FTOTWT,NUM,DUMA,DUMAX,DUMAY,DUMAZ,
     >                  DUMB,DUMBX,DUMBY,DUMBZ,AOX,GAOX,GAOY,GAOZ,G2AOXX
     *                 ,G2AOYY,G2AOZZ,G2AOXY,G2AOYZ,G2AOXZ,DA,DB,NEEDGR,
     *                  UROHF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL NEEDGR,UROHF
      PARAMETER(TWO=2.0D+00,ZERO=0.0D+00)
      PARAMETER (MXSH=5000, MXGTOT=20000)
      DIMENSION AOX(*),GAOX(*),GAOY(*),GAOZ(*),
     >     G2AOXX(*),G2AOYY(*),G2AOZZ(*),G2AOXY(*),G2AOYZ(*),G2AOXZ(*),
     >     DEDFT(3,*),DA(NUM,*),DB(NUM,*)
C
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
C***********************************************************************
C
C     COMPUTE DFT ENERGY GRADIENT (FUNCTIONAL PART ONLY)
C
C***********************************************************************
C
C     FOR A FUNCTIONAL FXC[RO,DELTA_RO] WHERE DELTA_RO IS GRADIENT OF
C     DENSITY TOTAL EXCHANGE-CORRELATION ENERGY E=INTEGRAL FXC DR
C     BY TAKING DERIVATIVE OVER SOME NUCLEAR COORDINATE R,
C     DE/DR =
C       INTEGRAL DF/D(RO)*D(RO)/DR   DF/D(DELTA_RO)*D(DELTA_RO)/DR DR
C     D(RO)/DR = - D(RO)/DR AND D(DELTA_RO)/DR = - D(DELTA_RO)/DR
C     (WHERE R IS ELECTRON COORDINATE)
C     DF/D(RO) AND DF/D(DELTA_RO) ARE STORED IN DUMA,DUMAX,DUMAY,DUMAZ
C     FOR ALPHA DENSITY AND IN DUMB, DUMBX, DUMBY, DUMBZ FOR BETA.
C     ALPHA + BETA DENSITY FORM ADDITIVE CONTRIBUTIONS TO THE GRADIENT.
C     NEXT,   D(RO)/DR(L) = 2 * SUM(I,J) D(I,J) DF(I)/DR(L) * F(J)
C     D(DELTA(K)_RO)/DR(L)= 2 * SUM(I,J) D(I,J)(DG[K](I)/DR[L]*F(J)
C                             +G[K](I)*G[L](J))
C     WHERE F(I) IS AO NUMBER I AND G[K](I) IS DF(I)/DR[K],
C     AND D(I,J) IS DENSITY, STORED IN DA AND DB, FOR ALPHA AND BETA.
C     F(I) IS STORED IN AOX, G[K](I) IN GAOX, GAOY, GAOZ,
C     DG[K](I)/DR[L] IN G2AOXX,G2AOYY,G2AOZZ,G2AOXY,G2AOYZ,G2AOXZ.
C
      DO 100 ISHELL=1,NSHELL
        IAT =KATOM(ISHELL)
        IF(IAT.EQ.NCNTR) GOTO 100
        MINI=KMIN (ISHELL)
        MAXI=KMAX (ISHELL)
        LOCI=KLOC (ISHELL)-MINI
        GRADX=ZERO
        GRADY=ZERO
        GRADZ=ZERO
        DO IBFN=MINI,MAXI
          ICHI=LOCI+IBFN
          DRA=ZERO
          DRB=ZERO
          IF(NEEDGR) THEN
            DRXA=ZERO
            DRYA=ZERO
            DRZA=ZERO
            DRXB=ZERO
            DRYB=ZERO
            DRZB=ZERO
            IF(UROHF) THEN
            DO JCHI=1,NUM
              DIJA=DA(ICHI,JCHI)
              DIJB=DB(ICHI,JCHI)
C
C         COMPUTE SUM(J) D(I,J) F(J)
C
              DRA= DRA+ DIJA*AOX(JCHI)
              DRB= DRB+ DIJB*AOX(JCHI)
C
C         COMPUTE SUM(J) D(I,J) G[K](J)
C
              DRXA=DRXA+DIJA*GAOX(JCHI)
              DRYA=DRYA+DIJA*GAOY(JCHI)
              DRZA=DRZA+DIJA*GAOZ(JCHI)
              DRXB=DRXB+DIJB*GAOX(JCHI)
              DRYB=DRYB+DIJB*GAOY(JCHI)
              DRZB=DRZB+DIJB*GAOZ(JCHI)
            ENDDO
            ELSE
            DO JCHI=1,NUM
              DIJA=DA(ICHI,JCHI)
C
C         COMPUTE SUM(J) D(I,J) F(J)
C
              DRA= DRA+ DIJA*AOX(JCHI)
C
C         COMPUTE SUM(J) D(I,J) G[K](J)
C
              DRXA=DRXA+DIJA*GAOX(JCHI)
              DRYA=DRYA+DIJA*GAOY(JCHI)
              DRZA=DRZA+DIJA*GAOZ(JCHI)
C
            ENDDO
            ENDIF
          ELSE
            IF(UROHF) THEN
              DO JCHI=1,NUM
                DRA=DRA+DA(ICHI,JCHI)*AOX(JCHI)
                DRB=DRB+DB(ICHI,JCHI)*AOX(JCHI)
              ENDDO
            ELSE
              DO JCHI=1,NUM
                DRA=DRA+DA(ICHI,JCHI)*AOX(JCHI)
              ENDDO
            ENDIF
          ENDIF
C
          GAOXI=GAOX(ICHI)
          GAOYI=GAOY(ICHI)
          GAOZI=GAOZ(ICHI)
C
C    COMPUTE -1/2*DF/D(RO)*D(RO)/DR=DF/D(RO)*SUM(I,J)D(I,J)DF(I)/DR*F(J)
C
          GRADX=GRADX+GAOXI*DRA*DUMA
          GRADY=GRADY+GAOYI*DRA*DUMA
          GRADZ=GRADZ+GAOZI*DRA*DUMA
          IF(UROHF) THEN
            GRADX=GRADX+GAOXI*DRB*DUMB
            GRADY=GRADY+GAOYI*DRB*DUMB
            GRADZ=GRADZ+GAOZI*DRB*DUMB
          ENDIF
          IF(NEEDGR) THEN
C
C         COMPUTE -1/2* SUM(K) DF/D(DELTA[K]_RO)*D(DELTA[K]_RO)/DR[L] =
C                 SUM(K)DF/D(DELTA[K]_RO)*SUM(I,J)D(I,J)
C    *                  (DG[K](I)/DR[L]*F(J)+G[K](I)*G[L](J))
C
            GRADX=GRADX+(G2AOXX(ICHI)*DRA + GAOXI*DRXA)*DUMAX +
     *                  (G2AOXY(ICHI)*DRA + GAOXI*DRYA)*DUMAY +
     *                  (G2AOXZ(ICHI)*DRA + GAOXI*DRZA)*DUMAZ
            GRADY=GRADY+(G2AOXY(ICHI)*DRA + GAOYI*DRXA)*DUMAX +
     *                  (G2AOYY(ICHI)*DRA + GAOYI*DRYA)*DUMAY +
     *                  (G2AOYZ(ICHI)*DRA + GAOYI*DRZA)*DUMAZ
            GRADZ=GRADZ+(G2AOXZ(ICHI)*DRA + GAOZI*DRXA)*DUMAX +
     *                  (G2AOYZ(ICHI)*DRA + GAOZI*DRYA)*DUMAY +
     *                  (G2AOZZ(ICHI)*DRA + GAOZI*DRZA)*DUMAZ
            IF(UROHF) THEN
            GRADX=GRADX+(G2AOXX(ICHI)*DRB + GAOXI*DRXB)*DUMBX +
     *                  (G2AOXY(ICHI)*DRB + GAOXI*DRYB)*DUMBY +
     *                  (G2AOXZ(ICHI)*DRB + GAOXI*DRZB)*DUMBZ
            GRADY=GRADY+(G2AOXY(ICHI)*DRB + GAOYI*DRXB)*DUMBX +
     *                  (G2AOYY(ICHI)*DRB + GAOYI*DRYB)*DUMBY +
     *                  (G2AOYZ(ICHI)*DRB + GAOYI*DRZB)*DUMBZ
            GRADZ=GRADZ+(G2AOXZ(ICHI)*DRB + GAOZI*DRXB)*DUMBX +
     *                  (G2AOYZ(ICHI)*DRB + GAOZI*DRYB)*DUMBY +
     *                  (G2AOZZ(ICHI)*DRB + GAOZI*DRZB)*DUMBZ
            ENDIF
          ENDIF
        ENDDO
        DEDFT(1,IAT)=DEDFT(1,IAT)-TWO*FTOTWT*GRADX
        DEDFT(2,IAT)=DEDFT(2,IAT)-TWO*FTOTWT*GRADY
        DEDFT(3,IAT)=DEDFT(3,IAT)-TWO*FTOTWT*GRADZ
C
C       USING TRANSLATIONAL INVARIANCE SUM(B) GRAD(B) EXC = 0
C       TO GET GRAD(A) E = - SUM(B.NE.A) GRAD(B) EXC
C       FOR A=NCNTR, SINCE THE ABOVE CODE IS NOT CORRECT IN THIS CASE.
C
        DEDFT(1,NCNTR)=DEDFT(1,NCNTR)+TWO*FTOTWT*GRADX
        DEDFT(2,NCNTR)=DEDFT(2,NCNTR)+TWO*FTOTWT*GRADY
        DEDFT(3,NCNTR)=DEDFT(3,NCNTR)+TWO*FTOTWT*GRADZ
  100 CONTINUE
C
      RETURN
      END
C*MODULE DFTGRD  *DECK DFTDER
      SUBROUTINE DFTDER
C***********************************************************************
C
C     ANALYTIC DFT GRADIENT (FOR THE GRID DFT).
C
C     --- WRITTEN BY ---
C     MUNEAKI KAMIYA, TAKAO TSUNEDA, SUSUMU YANAGISAWA, DMITRI FEDOROV
C
C                                                UNIVERSITY OF TOKYO
C
C***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      LOGICAL UROHF,OUT,GOPARR,DSKWRK,MASWRK
C
      COMMON /FMCOM / X(1)
      COMMON /GRAD  / DE(3,MXATM)
      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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA UHF,ROHF/8HUHF     ,8HROHF    /
      DATA DEBUG/8HDEBUG   /, DFTGRD/8HDERDFT  /,GRD1/8HGRD1    /
      DATA CHECK/8HCHECK   /
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3=L1*L1
C
      UROHF = SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
      OUT=(EXETYP.EQ.DFTGRD.OR.EXETYP.EQ.DEBUG.OR.EXETYP.EQ.GRD1).AND.
     *     MASWRK
C
      CALL GRDDFT(L2,NDFTEN)
C
      CALL VALFM(LOADFM)
C
      IDFTEN  = LOADFM  + 1
      IG2AXX  = IDFTEN  + NDFTEN
      IG2AYY  = IG2AXX  + NUM
      IG2AZZ  = IG2AYY  + NUM
      IG2AXY  = IG2AZZ  + NUM
      IG2AYZ  = IG2AXY  + NUM
      IG2AXZ  = IG2AYZ  + NUM
      ICOEFFA = IG2AXZ  + NUM
      ICOEFFB = ICOEFFA + NUM*NUM
      LDA     = ICOEFFB + NUM*NUM
      LDB     = LDA     + L2
      LDDA    = LDB     + L2
      LDDB    = LDDA    + L3
      LDEDFT  = LDDB    + L3
      IUVEC   = LDEDFT  + NAT*3
      IDWTINT = IUVEC   + NAT*3
      IDWTTOT = IDWTINT + 3*NAT*NAT
      IDTOTWT = IDWTTOT + 3*NAT
      LAST    = IDTOTWT + 3*NAT
C
      NEED = LAST - LOADFM - 1
C
C     ---- GET MEMORY ----
C
      CALL GETFM(NEED)
C     WRITE(6,*) 'WWWAHAHAA',NDFTEN,NUM,L2,L3,NAT
      IF(MASWRK) WRITE(IW,9000) NEED
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C     READ OR GENERATE DFT GRID. NOTE THAT THIS CAN BE A DIFFERENT GRID
C     COMPARED TO THE DFT ENERGY, BECAUSE THE SYMMETRY OF THE GRADIENT
C     COMPONENTS IS NOT THE SAME AS OF THE ENERGY, SO WE IN FACT TURN
C     OFF USING OCTANT SYMMETRY, BUT KEEP ATOMIC CENTRE SYMMETRY.
C     IN CASE OF C1 AND NON-ABELIAN GROUPS, OCTANT SYMMETRY IS NOT USED
C     SO THE GRID INFORMATIONS CAN BE READ.
C
      CALL DFTSET(X(IDFTEN),0,.TRUE.)
C
      CALL DAREAD(IDAF,IODA,X(ICOEFFA),NUM*NUM,15,0)
      CALL DAREAD(IDAF,IODA,X(LDA),L2,16,0)
      CALL EXPND(X(LDA),X(LDDA),L1,0)
      IF (UROHF) THEN
         CALL DAREAD(IDAF,IODA,X(ICOEFFB),NUM*NUM,19,0)
         CALL DAREAD(IDAF,IODA,X(LDB),L2,20,0)
         CALL EXPND(X(LDB),X(LDDB),L1,0)
      ENDIF
      CALL VCLR(X(LDEDFT),1,NAT*3)
      CALL DFTGRAD(X(IDFTEN),X(IG2AXX),X(IG2AYY),X(IG2AZZ),X(IG2AXY),
     *             X(IG2AYZ),X(IG2AXZ),X(ICOEFFA),X(ICOEFFB),X(LDDA),
     *             X(LDDB),X(LDEDFT),X(IUVEC),X(IDWTINT),X(IDWTTOT),
     *             X(IDTOTWT),L1,TOTELE,OUT)
C
      IF(OUT) WRITE(IW,9100)
      IF(OUT) CALL EGOUT(X(LDEDFT),NAT)
      CALL VADD(DE,1,X(LDEDFT),1,DE,1,3*NAT)
C
  800 CONTINUE
      CALL RETFM(NEED)
      CALL TIMIT(1)
      RETURN
C
 9000 FORMAT(/1X,'GRID DFT GRADIENT REQUIRES',I10,' WORDS OF MEMORY.')
 9100 FORMAT(/1X,'TOTAL DFT GRADIENT TERMS')
      END
