C 19 SEP 05 - MWS - ADD TRUE NUCLEAR CHARGE ARRAY TO INFOA COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL
C  7 SEP 04 - MWS - PAD COMMON BLOCK INTFIL
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C 16 JAN 04 - DGF - DFTINP: READ SWOFF AND GTHRE INPUTS
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 26 MAR 02 - MWS - DFTINP: ADD ARGUMENT TO PERMIT FORCING OF DFT
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C 26 JUN 01 - MWS - GRID WARNINGS, ALTER COMMON BLOCK WFNOPT
C 13 JUN 01 - TT,SY,MK,DGF - IMPLEMENT GRID-BASED DFT
C 26 OCT 00 - MWS - INTRODUCE MXAO PARAMETER
C 25 MAR 00 - KRG - NEW MODULE FOR GRID-FREE DFT COMPUTATIONS
C
C*MODULE DFT     *DECK DFTDRVR
      SUBROUTINE DFTDRVR(A,AFOCK,BFOCK,ADENSIN,BDENSIN,V)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL HISPIN,DIRSCF,FDIFF,PACK2E
C
      PARAMETER (MXSH=5000, MXATM=2000, MXGTOT=20000)
C
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      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 /OPTSCF/ DIRSCF,FDIFF
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DIMENSION A(*),AFOCK((NUM+NUM*NUM)/2),BFOCK((NUM+NUM*NUM)/2),
     * ADENSIN((NUM+NUM*NUM)/2),BDENSIN((NUM+NUM*NUM)/2),V(NUM,NUM)
C
      DATA    RHF  /8HRHF     /,ROHF  /8HROHF    /,UHF  /8HUHF     /
C
      L3 = NUM*NUM
      L4=NUM+NAUXFUN
      L6=L4*L4
C
C   THIS IS A GRID FREE DENSITY FUNCTIONAL CODE
C   FOR INFORMATION ON DFT IMPLEMENTATION SEE J CHEM PHYS 98, 5612(1993)
C   POPLE, GILL AND JOHNSON.
C   FOR INFORMATION ON GRID-FREE DFT SEE CHEM. PHYS. LETT. VOLUME 214
C   NUMBER 3,4 PAGES 397-401.  JAN ALMLOF AND YU CHENG ZHENG
C   A THEOCHEM ARTICLE BY THEM HAS BEEN PUBLISHED:
C   JOURNAL OF MOLECULAR STRUCTURE (THEOCHEM) VOLUME 388, PAGES 277-284
C   A JOURNAL OF CHEMICAL PHYSICS ARTICLE ON THIS CODE HAS PUBLISHED:
C   K.R. GLAESEMANN AND M.S. GORDON, "INVESTIGATION OF A GRID-FREE
C   DENSITY FUNCTIONAL THEORY (DFT) APPROACH", J. CHEM. PHYS., 108,
C   9959 (1998).
C   A JOURNAL OF CHEMICAL PHYSICS ARTICLES HAS BEEN ACCEPTED
C   K.R. GLAESEMANN AND M.S. GORDON, "EVALUATION OF GRADIENT
C   CORRECTIONS IN GRID-FREE DENSITY FUNCTIONAL THEORY", J. CHEM.
C   PHYS., IN PRESS.
C
C   FOR EACH FUNCTIONAL, BOTH AN ENERGY AND A FOCK TERM HAVE TO BE
C   CALCULATED.  THIS IS BECAUSE THE DOT PRODUCT OF THE FOCK MATRIX
C   AND THE DENSITY MATRIX IS NOT THE ENERGY.  THAT'S JUST A
C   CONVENIENCE IN HARTREE-FOCK THAT WE DON'T ENJOY.  EQ. 7 OF POPLE'S
C   PAPER IS VERY IMPORTANT TO REMEMBER.
C
C   EQUATIONS A5 AND A6 OF POPLES PAPER ARE WRONG
C
      IF (SCFTYP.EQ.RHF) THEN
         HISPIN=.FALSE.
      ELSE IF (SCFTYP.EQ.ROHF .OR. SCFTYP.EQ.UHF) THEN
         HISPIN=.TRUE.
      ELSE
         CALL ABRT
      END IF
C  THE DENSITY MATRIX
      IADENS=1
      IF(HISPIN) THEN
         IBDENS=IADENS+L3
      ELSE
         IBDENS=IADENS
      END IF
C  THE MATRIX REPRESENTATION OF THE DENSITY
      IALPHA=IBDENS+L3
      IF(HISPIN) THEN
        IBETA=IALPHA+L6
      ELSE
        IBETA=IALPHA
      END IF
C  LCAO COEFFICIENTS
      IVAUX=IBETA+L6
      IAUXATM=IVAUX+L6
      IAUXTYP=IAUXATM+NAUXSHL
      IAUXEXP=IAUXTYP+NAUXSHL
      IAUXCONT=IAUXEXP+NAUXSHL
      LOCAUX=IAUXCONT+NAUXSHL
      IMAPAUX=LOCAUX+NAUXSHL
C  THESE VARIABLES JUST KEEP GETTING USED
      ITEMP1=IMAPAUX+NAUXSHL*48
      ITEMP2=ITEMP1+L6
C
C  MEMORY FOR DFTRH3/4
C  (OVERWRITEN BY FUNCTIONALS, AFTER M[RHO] IS BUILT)
C
      IA3OV=ITEMP2+L6
      IND3OV=IA3OV+NINTMX
C  IF DIRECT SCF
      IF(DIRSCF .AND. (NSHELL*(NSHELL+1)/2) .GT. NINTMX) THEN
          IND3OV=IA3OV+NSHELL*(NSHELL+1)/2
      END IF
C  IF AUX FUNCTIONS
      IF(NAUXFUN .NE. 0 .AND. NSHELL*NSHELL .GT. NINTMX) THEN
          IND3OV=IA3OV+NSHELL*NSHELL
      END IF
C
      IA4OV=IA3OV
      IND4OV=IND3OV
C     NINTMX OR 15*15*15*15 FOR AUX FUNCTIONS
      IDD=IND3OV+15*15*15*15
      IDOP=IDD+NUM
      IDDVV=IDOP+NUM
      IDOPVV=IDDVV+NUM
C
      IAOAO=IDOPVV+NUM
      IAOAUX=IAOAO+(NSHELL+1)*NSHELL/2
      IAUXAUX=IAOAUX+NSHELL*NAUXSHL
      ITMPTMP=IAUXAUX+(NAUXSHL+1)*NAUXSHL/2
C
C     THESE TERMS OVERWRITE IA3OV ETC. MEMORY
C
      IF(DFTTYP(1) .GE. 2.0D+00) THEN
C  GRADIENT TERMS
        IDX  = ITEMP2 + L6
        IDY  = IDX  + L6
        IDZ  = IDY  + L6
C MEMORY FOR DIAGONALIATION ROUTINES
        IWRK=IDZ+L6
      ELSE
        IWRK=ITEMP2+L6
      END IF
      IIWRK=IWRK+L4*8
C MEMORY FOR STORING OUTPUT
      IRR     = IIWRK   + L4
      IRSA    = IRR     + L6
      IF(HISPIN) THEN
        IROPEN= IRSA    + L6
        IRSB  = IROPEN  + L6
      ELSE
        IROPEN= IRR
        IRSB  = IRSA
      END IF
C MEMORY FOR PASSING IN A COPY OF THE MATRIX REP OF THE DENSITY
      IRR2    = IRSB    + L6
      IRSA2   = IRR2    + L6
      IF(HISPIN) THEN
        IROPEN2=IRSA2   + L6
        IRSB2 = IROPEN2 + L6
      ELSE
        IROPEN2=IRR2
        IRSB2=IRSA2
      END IF
C MEMORY FOR FUNCTIONALS TO USE
      IWRKFN  = IRSB2   + L6
      ITEMP3  = IWRKFN  + L4*50
      ITEMP4  = ITEMP3  + L6
      ITEMP5  = ITEMP4  + L6
      ITEMP6  = ITEMP5  + L6
      ITEMP7  = ITEMP6  + L6
      ITEMP8  = ITEMP7  + L6
      ITEMP9  = ITEMP8  + L6
      ITEMP10 = ITEMP9  + L6
      ITEMP11 = ITEMP10 + L6
C
      IF (.NOT. HISPIN) THEN
         CALL CPYTSQ(ADENSIN,A(IADENS),NUM,1)
         CALL DSCAL(L3,0.5D+00,A(IADENS),1)
      ELSE
         CALL CPYTSQ(ADENSIN,A(IADENS),NUM,1)
         CALL CPYTSQ(BDENSIN,A(IBDENS),NUM,1)
      END IF
C
      IF(IDFT34.EQ.3) THEN
         CALL DFTRH3(A(IADENS),A(IBDENS),A(IA3OV),A(IND3OV),A(IDOP),
     *               A(IDD),A(IALPHA),A(IBETA),HISPIN,V,A(ITEMP1),
     *               A(ITEMP2),A(IDDVV),A(IDOPVV),NAUXFUN,NINTMX)
      ELSE IF(IDFT34.EQ.4) THEN
         CALL DFTRH4(A(IADENS),A(IBDENS),A(IA4OV),A(IALPHA),A(IBETA),
     *               HISPIN,A(ITEMP1),A(IND4OV),V,NAUXFUN,NINTMX)
      ELSE
         WRITE(IW,*)' IDFT NOT SET TO 3 OR 4 IN DFT MAIN DRIVER'
         CALL ABRT
      END IF
C
C     NOW AUXILIARY FUNCTIONS
C
      IF (NAUXFUN .EQ. 0) THEN
         CALL DCOPY(NUM*NUM,V,1,A(IVAUX),1)
       ELSE
C        COMBINE V AND VAUX
         CALL DFTRN4(V,A(IVAUX),NUM,NAUXFUN)
C        THROW VAUX ON THE END
         CALL DAREAD(IDAF,IODA,A(IVAUX+NUM*L4),NAUXFUN*L4,343,0)
C        THE R AND ROPEN MATRIX
         CALL DCOPY(NUM*NUM,A(IALPHA),1,A(ITEMP2),1)
         CALL DFTRN4(A(ITEMP2),A(IALPHA),NUM,NAUXFUN)
         IF (HISPIN) THEN
          CALL DCOPY(NUM*NUM,A(IBETA),1,A(ITEMP2),1)
          CALL DFTRN4(A(ITEMP2),A(IBETA),NUM,NAUXFUN)
         END IF
C
C        READ IN THE DESCRIPTION OF THE AUX BASIS SET
C        THIS SPILLS OVER IN MEMORY INTO THE RIGHT PLACES
C
         CALL DAREAD (IDAF,IODA,A(IAUXATM),NAUXSHL*53,349,0)
C
C        READ IN SCREENING INTEGRALS
C
         CALL DAREAD (IDAF,IODA,A(IAOAO),(NSHELL+1)*NSHELL/2,342,0)
         CALL DAREAD (IDAF,IODA,A(IAOAUX),NSHELL*NAUXSHL,353,0)
         CALL DAREAD (IDAF,IODA,A(IAUXAUX),(NAUXSHL+1)*NAUXSHL/2,354,0)
C
C        DO (IJK AUX) AND (IJ AUX AUX) INTEGRALS TO GET R AND ROPEN
C
         CALL DFTAX5A(A(ITEMP1),A(ITEMP2),A(IND3OV),NAUXFUN,NAUXSHL,
     *          A(IAUXATM),A(IAUXTYP),A(IAUXEXP),A(IAUXCONT),A(LOCAUX),
     *               A(IALPHA),A(IBETA),A(IADENS),A(IBDENS),A(IA3OV),
     *           HISPIN,A(IAOAO),A(IAOAUX))
         CALL DFTAX5B(A(ITEMP1),A(ITEMP2),A(IND3OV),NAUXFUN,NAUXSHL,
     *          A(IAUXATM),A(IAUXTYP),A(IAUXEXP),A(IAUXCONT),A(LOCAUX),
     *               A(IALPHA),A(IBETA),A(IA3OV),
     *            HISPIN,A(IAOAO),A(IAOAUX),A(IAUXAUX),A(IMAPAUX),
     *            A(ITMPTMP))
C
C        TRANSFORM R AND ROPEN TO THE MO BASIS SET
C
         CALL DFTTRN2(A(IALPHA),A(IVAUX),A(ITEMP2),L4)
         IF(HISPIN) CALL DFTTRN2(A(IBETA),A(IVAUX),A(ITEMP2),L4)
      END IF
C
      CALL VCLR(A(IRR),1,L6)
      CALL VCLR(A(IRSA),1,L6)
      IF(HISPIN) THEN
        CALL VCLR(A(IROPEN),1,L6)
        CALL VCLR(A(IRSB),1,L6)
      END IF
      IF(DFTTYP(1) .GE. 2.0D+00) THEN
         CALL DFTDX(A(IDX),A(IDY),A(IDZ),A(IVAUX),A(ITEMP1),NUM,NAUXFUN)
      END IF
C
C     NOW THAT WE HAVE THE MATRIX REPRESENTATION OF THE DENSITY
C     AND THE (I D/DX J) BOTH IN THE MO BASIS WE CAN DO SOME DFT
C
C     COPY THE MATRIX (ALPHA & BETA) TO TEMP MATRICES
C     (R2,RSA2,RSB2,ROPEN2) AND THEN DO DFT.  THEN ADD THE
C     INDIVIDUAL FUNCTIONAL TO TOTAL (R,RSA,RSB,ROPEN).
C
C     ALPHA EXCHANGE FUNCTIONALS
C
      IF(DFTTYP(2).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IALPHA),1,A(IRR2),1)
          CALL XALPHA(A(IRR2),A(IRSA2),A(ITEMP1),A(ITEMP2),A(ITEMP3),
     *                        A(IWRK),A(IIWRK),L4)
          CALL DAXPY(L6,DFTTYP(2),A(IRR2),1,A(IRR),1)
          CALL DAXPY(L6,DFTTYP(2),A(IRSA2),1,A(IRSA),1)
      END IF
C
C     DFTTYP(3) IS HARTREE-FOCK EXCHANGE
C
      IF(DFTTYP(4).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IALPHA),1,A(IRR2),1)
          CALL BECKE88(A(IRR2),A(IRSA2),A(IDX),A(IDY),A(IDZ),A(ITEMP1),
     *         A(ITEMP2),A(ITEMP3),A(ITEMP4),A(ITEMP5),A(ITEMP6),
     *         A(ITEMP7),A(ITEMP8),A(ITEMP9),A(ITEMP10),A(IWRK),
     *         A(IIWRK),L4)
          CALL DAXPY(L6,DFTTYP(4),A(IRR2),1,A(IRR),1)
          CALL DAXPY(L6,DFTTYP(4),A(IRSA2),1,A(IRSA),1)
      END IF
      IF(DFTTYP(5).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IALPHA),1,A(IRR2),1)
          CALL DEPRISTO(A(IRR2),A(IRSA2),A(IDX),A(IDY),A(IDZ),
     *         A(ITEMP1),
     *         A(ITEMP2),A(ITEMP3),A(ITEMP4),A(ITEMP5),A(ITEMP6),
     *         A(ITEMP7),
     *         A(ITEMP8),A(ITEMP9),A(ITEMP10),A(IWRK),A(IIWRK),
     *         L4)
          CALL DAXPY(L6,DFTTYP(5),A(IRR2),1,A(IRR),1)
          CALL DAXPY(L6,DFTTYP(5),A(IRSA2),1,A(IRSA),1)
      END IF
      IF(DFTTYP(6).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IALPHA),1,A(IRR2),1)
          CALL CAMA(A(IRR2),A(IRSA2),A(IDX),A(IDY),A(IDZ),A(ITEMP1),
     *         A(ITEMP2),A(ITEMP3),A(ITEMP4),A(ITEMP5),A(ITEMP6),
     *         A(ITEMP7),
     *         A(ITEMP8),A(ITEMP9),A(ITEMP10),A(IWRK),A(IIWRK),
     *         L4)
          CALL DAXPY(L6,DFTTYP(6),A(IRR2),1,A(IRR),1)
          CALL DAXPY(L6,DFTTYP(6),A(IRSA2),1,A(IRSA),1)
      END IF
      IF(DFTTYP(7).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IALPHA),1,A(IRR2),1)
          CALL CAMB(A(IRR2),A(IRSA2),A(IDX),A(IDY),A(IDZ),A(ITEMP1),
     *         A(ITEMP2),A(ITEMP3),A(ITEMP4),A(ITEMP5),A(ITEMP6),
     *         A(ITEMP7),A(ITEMP8),A(ITEMP9),A(ITEMP10),A(IWRK),
     *         A(IIWRK),A(ITEMP11),L4)
          CALL DAXPY(L6,DFTTYP(7),A(IRR2),1,A(IRR),1)
          CALL DAXPY(L6,DFTTYP(7),A(IRSA2),1,A(IRSA),1)
      END IF
C
C     BETA EXCHANGE FUNCTIONALS
C
      IF(HISPIN) THEN
      IF(DFTTYP(2).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IBETA),1,A(IROPEN2),1)
          CALL XALPHA(A(IROPEN2),A(IRSB2),A(ITEMP1),
     *            A(ITEMP2),A(ITEMP3),A(IWRK),A(IIWRK),L4)
          CALL DAXPY(L6,DFTTYP(2),A(IROPEN2),1,A(IROPEN),1)
          CALL DAXPY(L6,DFTTYP(2),A(IRSB2),1,A(IRSB),1)
      END IF
      IF(DFTTYP(4).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IBETA),1,A(IROPEN2),1)
          CALL BECKE88(A(IROPEN2),A(IRSB2),A(IDX),A(IDY),A(IDZ),
     *         A(ITEMP1),
     *         A(ITEMP2),A(ITEMP3),A(ITEMP4),A(ITEMP5),A(ITEMP6),
     *         A(ITEMP7),
     *         A(ITEMP8),A(ITEMP9),A(ITEMP10),A(IWRK),A(IIWRK),
     *         L4)
          CALL DAXPY(L6,DFTTYP(4),A(IROPEN2),1,A(IROPEN),1)
          CALL DAXPY(L6,DFTTYP(4),A(IRSB2),1,A(IRSB),1)
      END IF
      IF(DFTTYP(5).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IBETA),1,A(IROPEN2),1)
          CALL DEPRISTO(A(IROPEN2),A(IRSB2),A(IDX),A(IDY),A(IDZ),
     *         A(ITEMP1),
     *         A(ITEMP2),A(ITEMP3),A(ITEMP4),A(ITEMP5),A(ITEMP6),
     *         A(ITEMP7),
     *         A(ITEMP8),A(ITEMP9),A(ITEMP10),A(IWRK),A(IIWRK),
     *         L4)
          CALL DAXPY(L6,DFTTYP(5),A(IROPEN2),1,A(IROPEN),1)
          CALL DAXPY(L6,DFTTYP(5),A(IRSB2),1,A(IRSB),1)
      END IF
      IF(DFTTYP(6).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IBETA),1,A(IROPEN2),1)
          CALL CAMA(A(IROPEN2),A(IRSB2),A(IDX),A(IDY),A(IDZ),A(ITEMP1),
     *         A(ITEMP2),A(ITEMP3),A(ITEMP4),A(ITEMP5),A(ITEMP6),
     *         A(ITEMP7),
     *         A(ITEMP8),A(ITEMP9),A(ITEMP10),A(IWRK),A(IIWRK),
     *         L4)
          CALL DAXPY(L6,DFTTYP(6),A(IROPEN2),1,A(IROPEN),1)
          CALL DAXPY(L6,DFTTYP(6),A(IRSB2),1,A(IRSB),1)
      END IF
      IF(DFTTYP(7).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IBETA),1,A(IROPEN2),1)
          CALL CAMB(A(IROPEN2),A(IRSB2),A(IDX),A(IDY),A(IDZ),A(ITEMP1),
     *         A(ITEMP2),A(ITEMP3),A(ITEMP4),A(ITEMP5),A(ITEMP6),
     *         A(ITEMP7),A(ITEMP8),A(ITEMP9),A(ITEMP10),A(IWRK),
     *         A(IIWRK),A(ITEMP11),L4)
          CALL DAXPY(L6,DFTTYP(7),A(IROPEN2),1,A(IROPEN),1)
          CALL DAXPY(L6,DFTTYP(7),A(IRSB2),1,A(IRSB),1)
      END IF
      END IF
C
C    WIGNER LOCAL EXCHANGE-CORRELATION FUNCTIONALS
C
      IF(DFTTYP(8).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IALPHA),1,A(IRR2),1)
          IF(HISPIN) CALL DCOPY(L6,A(IBETA),1,A(IROPEN2),1)
          CALL WIGNER(A(IRR2),A(IRSA2),A(IROPEN2),A(IRSB2),
     *         A(ITEMP1),A(ITEMP2),A(ITEMP3),A(IWRK),A(IIWRK),HISPIN,
     *         L4)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(8),A(IROPEN2),1,A(IROPEN),1)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(8),A(IRSB2),1,A(IRSB),1)
          CALL DAXPY(L6,DFTTYP(8),A(IRR2),1,A(IRR),1)
          CALL DAXPY(L6,DFTTYP(8),A(IRSA2),1,A(IRSA),1)
      END IF
      IF(DFTTYP(9).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IALPHA),1,A(IRR2),1)
          IF(HISPIN) CALL DCOPY(L6,A(IBETA),1,A(IROPEN2),1)
          CALL WSDFT(A(IRR2),A(IRSA2),A(IROPEN2),A(IRSB2),
     *          A(ITEMP1),A(ITEMP2),A(ITEMP3),A(IWRK),A(IIWRK),HISPIN,
     *          L4)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(9),A(IROPEN2),1,A(IROPEN),1)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(9),A(IRSB2),1,A(IRSB),1)
          CALL DAXPY(L6,DFTTYP(9),A(IRR2),1,A(IRR),1)
          CALL DAXPY(L6,DFTTYP(9),A(IRSA2),1,A(IRSA),1)
      END IF
      IF(DFTTYP(10).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IALPHA),1,A(IRR2),1)
          IF(HISPIN) CALL DCOPY(L6,A(IBETA),1,A(IROPEN2),1)
          CALL WIGEXP(A(IRR2),A(IRSA2),A(IROPEN2),A(IRSB2),
     *          A(ITEMP1),A(ITEMP2),A(ITEMP3),A(IWRK),A(IIWRK),HISPIN,
     *          L4)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(10),A(IROPEN2),1,A(IROPEN),1)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(10),A(IRSB2),1,A(IRSB),1)
          CALL DAXPY(L6,DFTTYP(10),A(IRR2),1,A(IRR),1)
          CALL DAXPY(L6,DFTTYP(10),A(IRSA2),1,A(IRSA),1)
      END IF
C
C   LOCAL CORRELATION FUNCTIONALS
C
      IF(DFTTYP(11).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IALPHA),1,A(IRR2),1)
          IF(HISPIN) CALL DCOPY(L6,A(IBETA),1,A(IROPEN2),1)
          CALL VWN(A(IRR2),A(IRSA2),A(IROPEN2),A(IRSB2),A(ITEMP1),
     *               A(ITEMP2),A(ITEMP3),A(IWRK),A(IIWRK),HISPIN,
     *               A(ITEMP4),A(ITEMP5),A(ITEMP6),A(ITEMP7),A(ITEMP8),
     *               A(ITEMP9),L4)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(11),A(IROPEN2),1,A(IROPEN),1)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(11),A(IRSB2),1,A(IRSB),1)
          CALL DAXPY(L6,DFTTYP(11),A(IRR2),1,A(IRR),1)
          CALL DAXPY(L6,DFTTYP(11),A(IRSA2),1,A(IRSA),1)
      END IF
      IF(DFTTYP(12).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IALPHA),1,A(IRR2),1)
          IF(HISPIN) CALL DCOPY(L6,A(IBETA),1,A(IROPEN2),1)
          CALL PWLOC(A(IRR2),A(IRSA2),A(IROPEN2),A(IRSB2),A(ITEMP1),
     *               A(ITEMP2),A(ITEMP3),A(IWRK),A(IIWRK),HISPIN,
     *               A(ITEMP4),A(ITEMP5),A(ITEMP6),A(ITEMP7),A(ITEMP8),
     *               A(ITEMP9),L4)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(12),A(IROPEN2),1,A(IROPEN),1)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(12),A(IRSB2),1,A(IRSB),1)
          CALL DAXPY(L6,DFTTYP(12),A(IRR2),1,A(IRR),1)
          CALL DAXPY(L6,DFTTYP(12),A(IRSA2),1,A(IRSA),1)
      END IF
C
C    GRADIENT CORRECTED CORRELATION FUNCTIONALS
C
      IF(DFTTYP(13).NE.0.0D+00) THEN
          CALL DCOPY(L6,A(IALPHA),1,A(IRR2),1)
          IF(HISPIN) CALL DCOPY(L6,A(IBETA),1,A(IROPEN2),1)
          CALL LYP(A(IRR2),A(IRSA2),A(IROPEN2),A(IRSB2),
     *         A(IDX),A(IDY),A(IDZ),A(ITEMP1),A(ITEMP2),A(ITEMP3),
     *         A(ITEMP4),A(ITEMP5),A(ITEMP6),A(ITEMP7),A(ITEMP8),
     *         A(ITEMP9),A(ITEMP10),A(IWRK),A(IIWRK),A(ITEMP11),
     *         HISPIN,L4)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(13),A(IROPEN2),1,A(IROPEN),1)
          IF(HISPIN) CALL DAXPY(L6,DFTTYP(13),A(IRSB2),1,A(IRSB),1)
          CALL DAXPY(L6,DFTTYP(13),A(IRR2),1,A(IRR),1)
          CALL DAXPY(L6,DFTTYP(13),A(IRSA2),1,A(IRSA),1)
      END IF
C
C    NON-EXISTENT FUNCTIONALS
C
      IF(DFTTYP(14).NE.0.0D+00) THEN
          CALL ABRT
      END IF
      IF(DFTTYP(15).NE.0.0D+00) THEN
          CALL ABRT
      END IF
      IF(DFTTYP(16).NE.0.0D+00) THEN
          CALL ABRT
      END IF
      IF(DFTTYP(17).NE.0.0D+00) THEN
          CALL ABRT
      END IF
      IF(DFTTYP(18).NE.0.0D+00) THEN
          CALL ABRT
      END IF
      IF(DFTTYP(19).NE.0.0D+00) THEN
          CALL ABRT
      END IF
C
C     CALCULATE DFT ENERGY AND ADD DFT CONTRIBUTION TO THE FOCK MATRIX
C
      CALL DFTWRIT(A(IRR),A(IROPEN),A(IRSA),A(IRSB),
     *      A(ITEMP1),A(ITEMP2),A(IVAUX),AFOCK,BFOCK,A(ITEMP3),HISPIN,
     *      A(IADENS),A(IBDENS))
      RETURN
      END
C*MODULE DFT     *DECK DFTDX
      SUBROUTINE DFTDX(DX,DY,DZ,V,S,N,NAUXFUN)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DX(N+NAUXFUN,N+NAUXFUN),DY(N+NAUXFUN,N+NAUXFUN),
     *          DZ(N+NAUXFUN,N+NAUXFUN),V(N+NAUXFUN,N+NAUXFUN),S(*)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      L2=((N+NAUXFUN)+(N+NAUXFUN)*(N+NAUXFUN))/2
      IF(NAUXFUN .NE. 0) THEN
         IOFFSET=344
      ELSE
         IOFFSET=84
      END IF
C
C     PROCECESS THE DERIVATIVE MATRIX
C
      CALL DAREAD (IDAF,IODA,S,L2,0+IOFFSET,0)
      CALL EXPND  (S,DX,N+NAUXFUN,1)
      CALL DAREAD (IDAF,IODA,S,L2,1+IOFFSET,0)
      CALL EXPND  (S,DY,N+NAUXFUN,1)
      CALL DAREAD (IDAF,IODA,S,L2,2+IOFFSET,0)
      CALL EXPND  (S,DZ,N+NAUXFUN,1)
C
C   CONVERT TO MO BASIS
C
      CALL DFTTRN2(DX,V,S,N+NAUXFUN)
      CALL DFTTRN2(DY,V,S,N+NAUXFUN)
      CALL DFTTRN2(DZ,V,S,N+NAUXFUN)
      END
C*MODULE DFT     *DECK DFTINP
      SUBROUTINE DFTINP(FORCEDFT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION METHOD,MOROKM
C
      LOGICAL DSKWRK,MASWRK,GOPARR,THREETF
C
C      /DFTHF/,/DFTCF/ ARE INFORMATION FOR GRID CODE
C      /DFTPAR/ IS INFORMATION FOR GRID-FREE CODE
C      THE CALLING ARGUMENT SHOULD NORMALLY BE 'NONE', BUT MAY BE
C      ANY ONE OF THE LEGAL DFTTYP VALUES, IN ORDER TO FORCE DFT
C      CALCULATIONS (USEFUL FOR THE RUNTYP=SURFACE OPTION).
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 /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      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
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
C
      PARAMETER (ZERO=0.0D+0)
C
      PARAMETER (NNAM=16)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA DFT /8HDFT     /
      DATA QNAM/8HMETHOD  ,8HDFTTYP  ,8HTHREE   ,8HAUXFUN  ,8HNRAD    ,
     *          8HNTHE    ,8HNPHI    ,8HQOP     ,8HPFTYP   ,8HTHRESH  ,
     *          8HNRAD0   ,8HNTHE0   ,8HNPHI0   ,8HSWITCH  ,8HSWOFF   ,
     *          8HGTHRE   /
      DATA KQNAM /5,5,0,5,1,   1,1,3,5,3,  1,1,1,3,3,  3/
C
      DATA BECKE   /8HBECKE   /,DEPRIS /8HDEPRISTO/
      DATA HNONE   /8HNONE    /,CAMA   /8HCAMA    /,CAMB   /8HCAMB    /
      DATA WIGNER  /8HWIGNER  /,WS     /8HWS      /,WIGEXP /8HWIGEXP  /
      DATA VWN     /8HVWN     /,PWLOC  /8HPWLOC   /,HLYP   /8HLYP     /
      DATA BVWN    /8HBVWN    /,BLYP   /8HBLYP    /,B3LYP  /8HB3LYP   /
      DATA BPWLOC  /8HBPWLOC  /,XVWN   /8HXVWN    /,XPWLOC /8HXPWLOC  /
      DATA SLATER  /8HSLATER  /,SVWN   /8HSVWN    /,XALPHA /8HXALPHA  /
      DATA SPWLOC  /8HSPWLOC  /,HALF   /8HHALF    /
C
      DATA AUXBAS0 /8HAUX0    /,AUXBAS1 /8HAUX1    /
      DATA AUXBAS2 /8HAUX2    /,AUXBAS3 /8HAUX3    /
C
      DATA GRID/8HGRID    /
C
      DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF    /
      DATA HESS,MOROKM/8HHESSIAN ,8HMOROKUMA/
C
C     DEFAULT IS NO DFT, WHETHER GRID-FREE OR GRIDS
C
      DFTYPE=HNONE
      NAUXSHL=0
      NAUXFUN=0
      CALL VCLR(DFTTYP,1,20)
      NDFTFG=0
C
      METHOD = GRID
C
C     DEFAULT FOR GRID-FREE IS TO USE EXACT 4-CENTER METHOD INSTEAD OF
C     3-CENTER METHOD (NOT AVAILIBLE ANYWAY WITH AUXILIARY FUNCTIONS)
C
      AUXFUN = AUXBAS3
      THREETF=.FALSE.
C
C     DEFAULTS FOR GRID BASED CODE ARE
C
      PFTYP=HNONE
      NRAD     = 96
      NTHE     = 12
      NPHI     = 24
      QOP      = ZERO
      DFTTHR   = ZERO
      NRAD0    = 0
      NTHE0    = 0
      NPHI0    = 0
      SW0      = 3.0D-04
      SWOFF    = 5.0D-03
      DFTGTHR  = 1.0D+00
      IF(NFG.GT.0) DFTGTHR=10.0D+00
C
C     SWOFF: TURN OFF DFT UNTIL SCF CONVERGENCE FALLS BELOW SWOFF.
C            0.0 IS A SPECIAL VALUE TURNING THIS FEATURE OFF.
C            SWOFF IS SENSIBLE TO SET TO A VALUE LARGER THAN SW0.
C
      JRET = 0
      CALL NAMEIO(IR,JRET,DFT,NNAM,QNAM,KQNAM,
     *            METHOD,DFTYPE,THREETF,AUXFUN,NRAD,NTHE,NPHI,QOP,
     *            PFTYP,DFTTHR,NRAD0,NTHE0,NPHI0,SW0,SWOFF,DFTGTHR,
     *         0,0,0,    0,0,0,0,0,
     *     0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,
     *     0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0)
      IF(JRET .EQ. 2) THEN
         IF (MASWRK) WRITE (IW,*)' ERROR READING DFT GROUP'
         IF (MASWRK) WRITE (IW,*)' MAYBE YOU SPELLED DFTTYP AS DFTYP?'
         CALL ABRT
      END IF
C
      IF(FORCEDFT.NE.HNONE) DFTYPE=FORCEDFT
C
C        IF NOT DOING DFT, ITS TIME TO GET OUT OF HERE.
C
      IF(DFTYPE.EQ.HNONE) DFTTYP(3)=1.0D+00
      IF(DFTYPE.EQ.HNONE) RETURN
C                         ******
C
      NERR = 0
      IF(SCFTYP.NE.RHF .AND. SCFTYP.NE.UHF .AND. SCFTYP.NE.ROHF) THEN
        IF(MASWRK) WRITE(IW,9040) SCFTYP
        NERR=NERR+1
      END IF
      IF(RUNTYP.EQ.MOROKM) THEN
        IF(MASWRK) WRITE(IW,9045)
        NERR=NERR+1
      END IF
      IF(NERR.GT.0) CALL ABRT
C
C         GRID POINT CHECKS
C     THERE WILL BE ANOTHER CHECK OF NTHE AND NPHI LATER, SINCE
C     IF ANY SYMMETRY IS PRESENT SOME MULTIPLICITY IS REQUIRED.
C     (E.G. FOR D2H, NTHE MUST BE EVEN AND NPHI A MULTIPLE OF 4).
C
      IF(RUNTYP.EQ.HESS) THEN
         IF(NRAD0.EQ.0) NRAD0=NRAD
         IF(NPHI0.EQ.0) NPHI0=NPHI
         IF(NTHE0.EQ.0) NTHE0=NTHE
      END IF
      IF(NRAD0.EQ.0) NRAD0=MIN(MAX(NRAD/4,24),NRAD)
      IF(NTHE0.EQ.0) NTHE0=MIN(MAX(NTHE/3,8),NTHE)
      IF(NPHI0.EQ.0) NPHI0=MIN(MAX(NPHI/3,16),NPHI)
      IF(NRAD.LT.NRAD0.OR.NTHE.LT.NTHE0.OR.NPHI.LT.NPHI0) THEN
         IF(MASWRK) WRITE(IW,9060) NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0
         CALL ABRT
      ENDIF
      IF(NRAD.LT.1.OR.NTHE.LT.2.OR.NPHI.LT.1) THEN
         IF(MASWRK) WRITE(IW,9070) NRAD,NTHE,NPHI
         CALL ABRT
      ENDIF
      IF(NPHI.NE.2*NTHE  .OR.  NPHI0.NE.2*NTHE0) THEN
         IF(MASWRK) WRITE(IW,9080) NPHI,NTHE,NPHI0,NTHE0
      END IF
C
C      GRID DFT CALCULATION
C      (THIS CALL INCLUDES PRINTING ALL RELEVANT INPUT DATA)
C
      IF(METHOD.EQ.GRID) THEN
         CALL INPGDFT(DFTYPE,PFTYP)
         RETURN
      ENDIF
C
C            PROCESS AND PRINT OPTIONS FOR GRID-FREE DFT
C
      IF (MASWRK) THEN
         WRITE (IW,9020) DFTYPE,AUXFUN,THREETF
      END IF
C
      IF(DFTYPE.EQ.XALPHA) THEN
          IF(MASWRK) WRITE (IW,*)'    XALPHA-EXCHANGE (ALPHA=0.7)'
          DFTTYP(1)=1.0D+00
          DFTTYP(2)=1.05D+00
      ELSE IF(DFTYPE.EQ.SLATER) THEN
          IF(MASWRK) WRITE (IW,*)'    LSDA-EXCHANGE(SLATER)(ALPHA=2/3)'
          DFTTYP(1)=1.0D+00
          DFTTYP(2)=1.0D+00
      ELSE IF(DFTYPE.EQ.BECKE) THEN
          IF(MASWRK) WRITE (IW,*)'    BECKE88 EXCHANGE'
          DFTTYP(1)=2.0D+00
          DFTTYP(4)=1.0D+00
      ELSE IF(DFTYPE.EQ.DEPRIS) THEN
          IF(MASWRK) WRITE (IW,*)'    DEPRISTO-KRESS EXCHANGE'
          DFTTYP(1)=2.0D+00
          DFTTYP(5)=1.0D+00
      ELSE IF(DFTYPE.EQ.CAMA) THEN
          IF(MASWRK) WRITE (IW,*)'    CAMA EXCHANGE'
          DFTTYP(1)=2.0D+00
          DFTTYP(6)=1.0D+00
      ELSE IF(DFTYPE.EQ.CAMB) THEN
          IF(MASWRK) WRITE (IW,*)'    CAMB EXCHANGE'
          IF(MASWRK) WRITE (IW,*)'    CAMB CORRELATION'
          DFTTYP(1)=2.0D+00
          DFTTYP(7)=1.0D+00
      ELSE IF(DFTYPE.EQ.WIGNER) THEN
          IF(MASWRK) WRITE (IW,*)'    WIGNER EXCHANGE'
          IF(MASWRK) WRITE (IW,*)'    WIGNER CORRELATION'
          DFTTYP(1)=1.0D+00
          DFTTYP(8)=1.0D+00
      ELSE IF(DFTYPE.EQ.WS) THEN
          IF(MASWRK) WRITE (IW,*)'    WIGNER SCALED EXCHANGE'
          IF(MASWRK) WRITE (IW,*)'    WIGNER SCALED CORRELATION'
          DFTTYP(1)=1.0D+00
          DFTTYP(9)=1.0D+00
      ELSE IF(DFTYPE.EQ.WIGEXP) THEN
          IF(MASWRK) WRITE (IW,*)'    WIGNER EXPONENTIAL EXCHANGE'
          IF(MASWRK) WRITE (IW,*)'    WIGNER EXPONENTIAL CORRELATION'
          DFTTYP(1)=1.0D+00
          DFTTYP(10)=1.0D+00
      ELSE IF(DFTYPE.EQ.VWN) THEN
          IF(MASWRK) WRITE (IW,*)'    VWN CORRELATION'
          IF(MASWRK) WRITE (IW,*)'    HARTREE FOCK EXCHANGE'
          DFTTYP(1)=1.0D+00
          DFTTYP(11)=1.0D+00
          DFTTYP(3)=1.0D+00
      ELSE IF(DFTYPE.EQ.PWLOC) THEN
          IF(MASWRK) WRITE (IW,*)'    PW LOCAL CORRELATION'
          IF(MASWRK) WRITE (IW,*)'    HARTREE FOCK EXCHANGE'
          DFTTYP(1)=1.0D+00
          DFTTYP(12)=1.0D+00
          DFTTYP(3)=1.0D+00
      ELSE IF(DFTYPE.EQ.HLYP) THEN
          IF(MASWRK) WRITE (IW,*)
     *                       '    LYP CORRELATION (NOT SPIN POLARIZED)'
          IF(MASWRK) WRITE (IW,*)'    HARTREE FOCK EXCHANGE'
          DFTTYP(1)=2.0D+00
          DFTTYP(13)=1.0D+00
          DFTTYP(3)=1.0D+00
      ELSE IF(DFTYPE.EQ.BVWN) THEN
          IF(MASWRK) WRITE (IW,*)'    VWN CORRELATION'
          IF(MASWRK) WRITE (IW,*)'    BECKE88 EXCHANGE'
          DFTTYP(1)=2.0D+00
          DFTTYP(11)=1.0D+00
          DFTTYP(4)=1.0D+00
      ELSE IF(DFTYPE.EQ.BPWLOC) THEN
          IF(MASWRK) WRITE (IW,*)'    PW LOCAL CORRELATION'
          IF(MASWRK) WRITE (IW,*)'    BECKE88 EXCHANGE'
          DFTTYP(1)=2.0D+00
          DFTTYP(12)=1.0D+00
          DFTTYP(4)=1.0D+00
      ELSE IF(DFTYPE.EQ.XVWN) THEN
          IF(MASWRK) WRITE (IW,*)'    VWN LOCAL CORRELATION'
          IF(MASWRK) WRITE (IW,*)'    XALPHA-EXCHANGE (ALPHA=0.7)'
          DFTTYP(1)=1.0D+00
          DFTTYP(2)=1.05D+00
          DFTTYP(11)=1.0D+00
      ELSE IF(DFTYPE.EQ.SVWN) THEN
          IF(MASWRK) WRITE (IW,*)'    VWN LOCAL CORRELATION'
          IF(MASWRK) WRITE (IW,*)'    LSDA-EXCHANGE(SLATER)(ALPHA=2/3)'
          DFTTYP(1)=1.0D+00
          DFTTYP(2)=1.0D+00
          DFTTYP(11)=1.0D+00
      ELSE IF(DFTYPE.EQ.XPWLOC) THEN
          IF(MASWRK) WRITE (IW,*)'    PW LOCAL CORRELATION'
          IF(MASWRK) WRITE (IW,*)'    XALPHA-EXCHANGE (ALPHA=0.7)'
          DFTTYP(1)=1.0D+00
          DFTTYP(2)=1.05D+00
          DFTTYP(12)=1.0D+00
      ELSE IF(DFTYPE.EQ.SPWLOC) THEN
          IF(MASWRK) WRITE (IW,*)'    PW LOCAL CORRELATION'
          IF(MASWRK) WRITE (IW,*)'    LSDA-EXCHANGE(SLATER)(ALPHA=2/3)'
          DFTTYP(1)=1.0D+00
          DFTTYP(2)=1.0D+00
          DFTTYP(12)=1.0D+00
      ELSE IF(DFTYPE.EQ.BLYP) THEN
          IF(MASWRK) WRITE (IW,*)
     *                       '    LYP CORRELATION (NOT SPIN POLARIZED)'
          IF(MASWRK) WRITE (IW,*)'    BECKE88 EXCHANGE'
          DFTTYP(1)=2.0D+00
          DFTTYP(4)=1.0D+00
          DFTTYP(13)=1.0D+00
      ELSE IF(DFTYPE.EQ.HALF) THEN
          IF(MASWRK) WRITE (IW,*)'    HALF BECKE88 EXCHANGE'
          IF(MASWRK) WRITE (IW,*)'    HALF HARTREE FOCK EXCHANGE'
          DFTTYP(1)=2.0D+00
          DFTTYP(4)=0.5D+00
          DFTTYP(3)=0.5D+00
      ELSE IF(DFTYPE.EQ.B3LYP) THEN
          IF(MASWRK) WRITE (IW,*)'    B3LYP FUNCTIONAL'
          IF(MASWRK) WRITE (IW,*)'    USING VWN5 FOR LOCAL CORRELATION'
          IF(MASWRK) WRITE (IW,*)
     *                     '    LYP CORRELATION IS NOT SPIN POLARIZED'
          DFTTYP(1)=2.0D+00
          DFTTYP(2)=0.08D+00
          DFTTYP(3)=0.2D+00
          DFTTYP(4)=0.72D+00
          DFTTYP(11)=0.19D+00
          DFTTYP(13)=0.81D+00
      ELSE
          IF(MASWRK) WRITE (IW,*)'    INVALID DFT FUNCTIONAL'
          CALL ABRT
      END IF
C
      IF(THREETF) THEN
          IDFT34=3
          IF(MASWRK) WRITE (IW,*)'    WILL USE 3 CENTER INTEGRALS'
      ELSE
          IDFT34=4
          IF(MASWRK) WRITE (IW,*)'    WILL USE 4 CENTER INTEGRALS'
      END IF
C
      IF (AUXFUN .EQ. HNONE .OR. AUXFUN .EQ. AUXBAS0) THEN
          IF(MASWRK) WRITE (IW,*)'    USING NO AUXILIARY BASIS SET!!'
          IF(MASWRK) WRITE (IW,*)'    I SMELL TROUBLE BREWING'
      ELSE IF (AUXFUN .EQ. AUXBAS1) THEN
          DFTTYP(20)=1.0D+00
          IF(MASWRK) WRITE (IW,*)'    USING AUX1 AUXILIARY BASIS SET'
      ELSE IF (AUXFUN .EQ. AUXBAS2) THEN
          DFTTYP(20)=2.0D+00
          IF(MASWRK) WRITE (IW,*)'    USING AUX2 AUXILIARY BASIS SET'
      ELSE IF (AUXFUN .EQ. AUXBAS3) THEN
          DFTTYP(20)=3.0D+00
          IF(MASWRK) WRITE (IW,*)'    USING AUX3 AUXILIARY BASIS SET'
      ELSE
          IF(MASWRK) WRITE (IW,*)'    ERROR IN AUXILIARY BASIS INPUT'
          CALL ABRT
      END IF
C
      IF (AUXFUN .NE. HNONE .AND. IDFT34 .EQ. 3) THEN
          IF(MASWRK) WRITE (IW,*)'    THREE=.T. MUST HAVE AUXFUN=NONE'
          CALL ABRT
      END IF
C
C      THE DFTTYP ARRAY CONSISTS OF
C      DFTTYP(1)=0 NO DFT
C      DFTTYP(1)=1 LOCAL DFT
C      DFTTYP(1)=2 GRADIENT CORRECTED DFT
C      DFTTYP(2-19)= FRACTIONAL WEIGHTS OF FUNCTIONALS
C      2  X-ALPHA WITH A=2/3
C      3  HARTREE-FOCK EXACT EXCHANGE
C      4  BECKE 88 EXCHANGE
C      5  DE PRISTO-KRESS EXCHANGE
C      6  CAMA EXCHANGE
C      7  CAMB EXCHANGE (AND CORRELATION)
C      8  WIGNER EXCHANGE-CORRELATION
C      9  WIGNER SCALED FORM
C     10  WIGNER EXPONENTIAL FORM
C     11  VWN5 LOCAL CORRELATION
C     12  PERDEW-WANG LOCAL CORRELATION
C     13  LYP GRADIENT CORRECTED CORRELATION (NOT SPIN POLARIZED)
C     14  REMAINING ARE UNUSED
C     15
C     16
C     17
C     18
C     19
C     20  TYPE OF AUXILIARY BASIS SET USED
C
      RETURN
C
C            PROCESS AND PRINT OPTIONS FOR GRID-FREE DFT
C     WRITTEN BY TAKAO TSUNEDA, UNIVERSITY OF TOKYO
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
 9020 FORMAT(/5X,'GRID-FREE DFT OPTIONS'/5X,15(1H-)/
     *  5X,7HDFTTYP=,A8,5X,7HAUXFUN=,A8,5X,7HTHREE =,L8)
 9040 FORMAT(1X,'ERROR, DFT IS AVAILABLE FOR SCFTYP=RHF, UHF, OR ROHF'/
     *       1X,'BUT YOU CHOSE SCFTYP=',A8)
 9045 FORMAT(1X,'DFT MAY NOT BE USED WITH RUNTYP=MOROKUMA')
 9060 FORMAT(1X,'NRAD,NPHI,NTHE SHOULD DEFINE A FINER GRID THAN ',
     *          'NRAD0,NPHI0,NTHE0',/1X,6I4)
 9070 FORMAT(1X,'NRAD,NPHI,NTHE ARE INVALID.',3I4)
 9080 FORMAT(//10X,'* * * * WARNING * * * *'/
     *       1X,'DFT ANGULAR GRID DOES NOT HAVE PHI TWICE THETA.'/
     *       1X,'THIS IS NOT VERY SENSIBLE IN SPHERICAL COORDINATES!'/
     *       1X,'NPHI,NTHE=',2I4,' NPHI0,NTHE0=',2I4//)
      END
C*MODULE DFT     *DECK DFTMEM
      SUBROUTINE DFTMEM(IDFTMEM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXSH=5000,MXATM=2000, MXGTOT=20000)
C
      LOGICAL DIRSCF,FDIFF,HISPIN,PACK2E
C
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      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 /OPTSCF/ DIRSCF,FDIFF
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA    RHF  /8HRHF     /
C
      NUMMY=NUM+NAUXFUN
      NEED1=0
      NEED2=0
      MAX=3
      L3 = NUMMY*NUMMY
      IF(SCFTYP .EQ. RHF) THEN
        HISPIN=.FALSE.
      ELSE
        HISPIN=.TRUE.
      END IF
      IF(HISPIN)NEED1=NEED1+2*L3
      NEED1=NEED1+NINTMX+15*15*15*15+4*NUMMY+4*L3
      IF(DIRSCF .AND. (NSHELL*(NSHELL+1)/2).GT.15*15*15*15 .AND.
     *                      NAUXFUN.EQ.0) THEN
          NEED1=NEED1+NSHELL*(NSHELL+1)/2-15*15*15*15
      ELSE IF (NAUXFUN .NE. 0 .AND. NSHELL*NSHELL .GT. 15*15*15*15) THEN
          NEED1=NEED1+NSHELL*NSHELL-15*15*15*15
      END IF
      IF(NAUXFUN .NE. 0) THEN
C         SCREENING INTS
          NEED1=NEED1+((NSHELL+NAUXSHL+1)*(NSHELL+NAUXSHL))/2
C         SYMMETRY MAPPING
          NEED1=NEED1+48*NAUXSHL+(NAUXSHL+1)*NAUXSHL/2
      END IF
C
      NEED2=NEED2+4*L3
      IF(HISPIN)NEED2=NEED2+2*L3
      IF(DFTTYP(1) .GE. 2.0D+00) NEED2=NEED2+3*L3
      NEED2=NEED2+9*NUMMY
      NEED2=NEED2+4*L3
      IF(HISPIN)NEED2=NEED2+4*L3
      NEED2=NEED2+50*NUMMY
      NEED2=NEED2+L3
C
      IF(DFTTYP(4).NE.0.0D+00 .AND. MAX .LT. 10) MAX=10
      IF(DFTTYP(5).NE.0.0D+00 .AND. MAX .LT. 10) MAX=10
      IF(DFTTYP(6).NE.0.0D+00 .AND. MAX .LT. 10) MAX=10
      IF(DFTTYP(7).NE.0.0D+00 .AND. MAX .LT. 11) MAX=11
      IF(DFTTYP(11).NE.0.0D+00 .AND. MAX .LT. 9) MAX=9
      IF(DFTTYP(12).NE.0.0D+00 .AND. MAX .LT. 9) MAX=9
      IF(DFTTYP(13).NE.0.0D+00 .AND. MAX .LT. 11) MAX=11
      NEED2=NEED2+(MAX-3)*L3
      IF (NEED2 .GT. NEED1) THEN
        IDFTMEM=NEED2
      ELSE
        IDFTMEM=NEED1
      END IF
C
C     ADD IN MEMORY FOR V WITH AUX FUNCTIONS
C
      IDFTMEM=IDFTMEM+L3+NAUXSHL*5
      END
C*MODULE DFT     *DECK DFTRH3
      SUBROUTINE DFTRH3(ADENS,BDENS,A3OV,IND3OV,DOP,DD,R,ROPEN,
     *                  HISPIN,V,S,VV,DDVV,DOPVV,NAUXFUN,NINTMX)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL HISPIN,SVDSKW,GOPARR,DSKWRK,MASWRK,DIRSCF,FDIFF,DONE
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IJPAIR/ IJINDEX(MXAO)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
      COMMON /OPTSCF/ DIRSCF,FDIFF
C
      DIMENSION ADENS(NUM,NUM),BDENS(NUM,NUM),A3OV(NINTMX),IND3OV(*),
     *          DOP(NUM),DD(NUM),R(NUM,NUM),ROPEN(NUM,NUM),
     *          V(NUM,NUM),S(NUM,NUM),VV(NUM,NUM),DDVV(NUM),DOPVV(NUM)
C
      SVDSKW = DSKWRK
      DSKWRK = .TRUE.
      IF(.NOT. DIRSCF) CALL SEQREW(21)
      DONE=.FALSE.
      L3=NUM*NUM
      IF (HISPIN) THEN
        CALL VCLR(ROPEN,1,L3)
        CALL VCLR(DOP,1,NUM)
      END IF
      CALL VCLR(R,1,L3)
      CALL VCLR(DD,1,NUM)
C
C     MAKE THE MATRIX REPRESENTATION OF THE TOTAL DENSITY
C     R(I,J) = <I|(DENSITY)|J>
C
C
      IF (.NOT. HISPIN) THEN
C       HALF THE DIAGONAL ELEMENTS
C       SO (IJKK) ISN'T COUNTED TWICE BUT (IJKL) IS
        DO 250 I=1,NUM
           ADENS(I,I)=ADENS(I,I)*0.5D+00
 250    CONTINUE
C
C
        IF(.NOT. DIRSCF) THEN
C
C
 555    CALL PREAD(21,A3OV,IND3OV,ILOOP,NINTMX)
        IF(ILOOP .LE. 0) THEN
           DONE=.TRUE.
           ILOOP = - ILOOP
        END IF
C  F77 MEANS THAT THIS LOOP WILL NOT EXECUTE IF ILOOP=ZERO
        DO 565 I=1,ILOOP
C
C   USE LABSIZ TO FIND POSITION (NUM AND BIT SIZE DEPENDANT)
C
                       NPACK = I
                       IF (LABSIZ .EQ. 2) THEN
*I32                     LABEL = IND3OV( 2*NPACK - 1 )
*I32                     IPACK = ISHFT( LABEL, -16 )
*I32                     JPACK = IAND(  LABEL, 65535 )
*I32                     LABEL = IND3OV( 2*NPACK     )
*I32                     KPACK = ISHFT( LABEL, -16 )
*I64                     LABEL = IND3OV(NPACK)
*I64                     IPACK = ISHFT( LABEL, -48 )
*I64                     JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                     KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
                       ELSE IF (LABSIZ .EQ. 1) THEN
*I32                     LABEL = IND3OV(NPACK)
*I32                     IPACK = ISHFT( LABEL, -24 )
*I32                     JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32                     KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                     IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                       LABEL = IND3OV( NPACK/2 )
*I64                       IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64                       JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                     ELSE
*I64                       LABEL = IND3OV( (NPACK/2)+1 )
*I64                       IPACK = ISHFT( LABEL, -56 )
*I64                       JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64                     END IF
                       END IF
                       K = IPACK
                       L = JPACK
                       N = KPACK
C
         A3OVKLN=A3OV(I)
         IF(K .NE. L .AND. L .NE. N) THEN
C            (KLN)
             DD(K)=DD(K)+A3OVKLN*ADENS(L,N)
             DD(L)=DD(L)+A3OVKLN*ADENS(K,N)
             DD(N)=DD(N)+A3OVKLN*ADENS(K,L)
         ELSEIF(K .NE. L) THEN
C            (KLL)
             DD(K)=DD(K)+A3OVKLN*ADENS(L,L)
             DD(L)=DD(L)+A3OVKLN*ADENS(K,L)
         ELSEIF(L .NE. N) THEN
C            (KKN)
             DD(K)=DD(K)+A3OVKLN*ADENS(K,N)
             DD(N)=DD(N)+A3OVKLN*ADENS(K,K)
         ELSE
C            (KKK)
             DD(K)=DD(K)+A3OVKLN*ADENS(K,K)
         END IF
C
 565      CONTINUE
      IF( .NOT. DONE) GO TO 555
C
C
      ELSE
C
C
      CALL THRDIR(IND3OV,HISPIN,.TRUE.,R,ROPEN,ADENS,BDENS,DD,DDVV,
     *            DOPVV,DOP)
C
C
      END IF
C
C
      IF(GOPARR) CALL DDI_GSUMF(1038,DD,NUM)
C
C     CALCULATE AO-->MO FACTOR
C
      DO 200 I=1,NUM
        DO 201 J=1,I
          DUM=0.0D+00
          DO 202 K=1,NUM
            DUM=DUM+V(J,K)*V(I,K)
 202      CONTINUE
          VV(I,J)=DUM
          VV(J,I)=DUM
 201    CONTINUE
 200  CONTINUE
C
C     USE THE AO TO MO FACTOR
C
      DO 210 I=1,NUM
        DUM=0.0D+00
        DO 211 J=1,NUM
          DUM=DUM+DD(J)*VV(J,I)
 211    CONTINUE
        DDVV(I)=DUM
 210  CONTINUE
C
C  WE CALCULATE MPQ[RHO]
C
        DONE=.FALSE.
C
C
        IF(.NOT. DIRSCF) THEN
C
C
        CALL SEQREW(21)
 550    CALL PREAD(21,A3OV,IND3OV,ILOOP,NINTMX)
        IF(ILOOP .LE. 0) THEN
           DONE=.TRUE.
           ILOOP = - ILOOP
        END IF
C  F77 MEANS THAT THIS LOOP WILL NOT EXECUTE IF ILOOP=ZERO
        DO 590 I=1,ILOOP
C
C   USE LABSIZ TO FIND POSITION (NUM AND BIT SIZE DEPENDANT)
C
                       NPACK = I
                       IF (LABSIZ .EQ. 2) THEN
*I32                     LABEL = IND3OV( 2*NPACK - 1 )
*I32                     IPACK = ISHFT( LABEL, -16 )
*I32                     JPACK = IAND(  LABEL, 65535 )
*I32                     LABEL = IND3OV( 2*NPACK     )
*I32                     KPACK = ISHFT( LABEL, -16 )
*I64                     LABEL = IND3OV(NPACK)
*I64                     IPACK = ISHFT( LABEL, -48 )
*I64                     JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                     KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
                       ELSE IF (LABSIZ .EQ. 1) THEN
*I32                     LABEL = IND3OV(NPACK)
*I32                     IPACK = ISHFT( LABEL, -24 )
*I32                     JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32                     KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                     IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                       LABEL = IND3OV( NPACK/2 )
*I64                       IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64                       JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                     ELSE
*I64                       LABEL = IND3OV( (NPACK/2)+1 )
*I64                       IPACK = ISHFT( LABEL, -56 )
*I64                       JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64                     END IF
                       END IF
                       K = IPACK
                       L = JPACK
                       N = KPACK
C
         A3OVKLN=A3OV(I)
         IF(K .NE. L .AND. L .NE. N) THEN
C            (KLN)
             R(K,L)=R(K,L)+A3OVKLN*DDVV(N)
             R(K,N)=R(K,N)+A3OVKLN*DDVV(L)
             R(L,N)=R(L,N)+A3OVKLN*DDVV(K)
         ELSEIF(K .NE. L) THEN
C            (KLL)
             R(K,L)=R(K,L)+A3OVKLN*DDVV(L)
             R(L,L)=R(L,L)+A3OVKLN*DDVV(K)
         ELSEIF(L .NE. N) THEN
C            (KKN)
             R(K,K)=R(K,K)+A3OVKLN*DDVV(N)
             R(K,N)=R(K,N)+A3OVKLN*DDVV(K)
         ELSE
C            (KKK)
             R(K,K)=R(K,K)+A3OVKLN*DDVV(K)
         END IF
 590    CONTINUE
        IF(.NOT. DONE) GO TO 550
C
C
      ELSE
C
C
      CALL THRDIR(IND3OV,HISPIN,.FALSE.,R,ROPEN,ADENS,BDENS,DD,DDVV,
     *            DOPVV,DOP)
C
C
      END IF
C
C       FILL OUT THE TRIANGLE AND SYMMETRIZE
C
        CALL DSCAL(L3,2.0D+00,R,1)
        IF(GOPARR) CALL DDI_GSUMF(1030,R,L3)
C
        CALL CPYSQT(R,S,NUM,1)
        CALL SYMH(S,R,IJINDEX)
        CALL EXPND(S,R,NUM,0)
        IF(NAUXFUN .EQ. 0) CALL DFTTRN2(R,V,S,NUM)
C
C       RESTORE THE DENSITY MATRIX
C
        DO 251 I=1,NUM
           ADENS(I,I)=ADENS(I,I)*2.0D+00
 251    CONTINUE
      ELSE
C
C   OPEN SHELL CASE
C
C
C       HALF THE DIAGONAL ELEMENTS
C       SO (IJKK) ISN'T COUNTED TWICE BUT (IJKL) IS
        DO 252 I=1,NUM
           ADENS(I,I)=ADENS(I,I)*0.5D+00
           BDENS(I,I)=BDENS(I,I)*0.5D+00
 252    CONTINUE
C
C
        IF(.NOT. DIRSCF) THEN
C
C
 655    CALL PREAD(21,A3OV,IND3OV,ILOOP,NINTMX)
        IF(ILOOP .LE. 0) THEN
           DONE=.TRUE.
           ILOOP = - ILOOP
        END IF
C  F77 MEANS THAT THIS LOOP WILL NOT EXECUTE IF ILOOP=ZERO
        DO 665 I=1,ILOOP
C
C   USE LABSIZ TO FIND POSITION (NUM AND BIT SIZE DEPENDANT)
C
                       NPACK = I
                       IF (LABSIZ .EQ. 2) THEN
*I32                     LABEL = IND3OV( 2*NPACK - 1 )
*I32                     IPACK = ISHFT( LABEL, -16 )
*I32                     JPACK = IAND(  LABEL, 65535 )
*I32                     LABEL = IND3OV( 2*NPACK     )
*I32                     KPACK = ISHFT( LABEL, -16 )
*I64                     LABEL = IND3OV(NPACK)
*I64                     IPACK = ISHFT( LABEL, -48 )
*I64                     JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                     KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
                       ELSE IF (LABSIZ .EQ. 1) THEN
*I32                     LABEL = IND3OV(NPACK)
*I32                     IPACK = ISHFT( LABEL, -24 )
*I32                     JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32                     KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                     IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                       LABEL = IND3OV( NPACK/2 )
*I64                       IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64                       JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                     ELSE
*I64                       LABEL = IND3OV( (NPACK/2)+1 )
*I64                       IPACK = ISHFT( LABEL, -56 )
*I64                       JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64                     END IF
                       END IF
                       K = IPACK
                       L = JPACK
                       N = KPACK
C
         A3OVKLN=A3OV(I)
         IF(K .NE. L .AND. L .NE. N) THEN
C            (KLN)
             DD(K)=DD(K)+A3OVKLN*ADENS(L,N)
             DD(L)=DD(L)+A3OVKLN*ADENS(K,N)
             DD(N)=DD(N)+A3OVKLN*ADENS(K,L)
             DOP(K)=DOP(K)+A3OVKLN*BDENS(L,N)
             DOP(L)=DOP(L)+A3OVKLN*BDENS(K,N)
             DOP(N)=DOP(N)+A3OVKLN*BDENS(K,L)
         ELSEIF(K .NE. L) THEN
C            (KLL)
             DD(K)=DD(K)+A3OVKLN*ADENS(L,L)
             DD(L)=DD(L)+A3OVKLN*ADENS(K,L)
             DOP(K)=DOP(K)+A3OVKLN*BDENS(L,L)
             DOP(L)=DOP(L)+A3OVKLN*BDENS(K,L)
         ELSEIF(L .NE. N) THEN
C            (KKN)
             DD(K)=DD(K)+A3OVKLN*ADENS(K,N)
             DD(N)=DD(N)+A3OVKLN*ADENS(K,K)
             DOP(K)=DOP(K)+A3OVKLN*BDENS(K,N)
             DOP(N)=DOP(N)+A3OVKLN*BDENS(K,K)
         ELSE
C            (KKK)
             DD(K)=DD(K)+A3OVKLN*ADENS(K,K)
             DOP(K)=DOP(K)+A3OVKLN*BDENS(K,K)
         END IF
C
 665      CONTINUE
      IF( .NOT. DONE) GO TO 655
C
C
      ELSE
C
C
      CALL THRDIR(IND3OV,HISPIN,.TRUE.,R,ROPEN,ADENS,BDENS,DD,DDVV,
     *            DOPVV,DOP)
C
C
      END IF
C
C
      IF(GOPARR) CALL DDI_GSUMF(1036,DD,NUM)
      IF(GOPARR) CALL DDI_GSUMF(1037,DOP,NUM)
C
C     CALCULATE AO-->MO FACTOR
C
      DO 300 I=1,NUM
        DO 301 J=1,I
          DUM=0.0D+00
          DO 302 K=1,NUM
            DUM=DUM+V(J,K)*V(I,K)
 302      CONTINUE
          VV(I,J)=DUM
          VV(J,I)=DUM
 301    CONTINUE
 300  CONTINUE
C
C     USE THE AO TO MO FACTOR
C
      DO 310 I=1,NUM
        DUM=0.0D+00
        DUM2=0.0D+00
        DO 311 J=1,NUM
          DUM=DUM+DD(J)*VV(J,I)
          DUM2=DUM2+DOP(J)*VV(J,I)
 311    CONTINUE
        DDVV(I)=DUM
        DOPVV(I)=DUM2
 310  CONTINUE
C
C  WE CALCULATE MPQ[RHO]
C
C
C
        IF(.NOT. DIRSCF) THEN
C
C
        DONE=.FALSE.
        CALL SEQREW(21)
 650    CALL PREAD(21,A3OV,IND3OV,ILOOP,NINTMX)
        IF(ILOOP .LE. 0) THEN
           DONE=.TRUE.
           ILOOP = - ILOOP
        END IF
C  F77 MEANS THAT THIS LOOP WILL NOT EXECUTE IF ILOOP=ZERO
        DO 680 I=1,ILOOP
C
C   USE LABSIZ TO FIND POSITION (NUM AND BIT SIZE DEPENDANT)
C
                       NPACK = I
                       IF (LABSIZ .EQ. 2) THEN
*I32                     LABEL = IND3OV( 2*NPACK - 1 )
*I32                     IPACK = ISHFT( LABEL, -16 )
*I32                     JPACK = IAND(  LABEL, 65535 )
*I32                     LABEL = IND3OV( 2*NPACK     )
*I32                     KPACK = ISHFT( LABEL, -16 )
*I64                     LABEL = IND3OV(NPACK)
*I64                     IPACK = ISHFT( LABEL, -48 )
*I64                     JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                     KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
                       ELSE IF (LABSIZ .EQ. 1) THEN
*I32                     LABEL = IND3OV(NPACK)
*I32                     IPACK = ISHFT( LABEL, -24 )
*I32                     JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32                     KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                     IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                       LABEL = IND3OV( NPACK/2 )
*I64                       IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64                       JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                     ELSE
*I64                       LABEL = IND3OV( (NPACK/2)+1 )
*I64                       IPACK = ISHFT( LABEL, -56 )
*I64                       JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64                     END IF
                       END IF
                       K = IPACK
                       L = JPACK
                       N = KPACK
C
         A3OVKLN=A3OV(I)
         IF(K .NE. L .AND. L .NE. N) THEN
C            (KLN)
             R(K,L)=R(K,L)+A3OVKLN*DDVV(N)
             R(K,N)=R(K,N)+A3OVKLN*DDVV(L)
             R(L,N)=R(L,N)+A3OVKLN*DDVV(K)
             ROPEN(K,L)=ROPEN(K,L)+A3OVKLN*DOPVV(N)
             ROPEN(K,N)=ROPEN(K,N)+A3OVKLN*DOPVV(L)
             ROPEN(L,N)=ROPEN(L,N)+A3OVKLN*DOPVV(K)
         ELSEIF(K .NE. L) THEN
C            (KLL)
             R(K,L)=R(K,L)+A3OVKLN*DDVV(L)
             R(L,L)=R(L,L)+A3OVKLN*DDVV(K)
             ROPEN(K,L)=ROPEN(K,L)+A3OVKLN*DOPVV(L)
             ROPEN(L,L)=ROPEN(L,L)+A3OVKLN*DOPVV(K)
         ELSEIF(L .NE. N) THEN
C            (KKN)
             R(K,K)=R(K,K)+A3OVKLN*DDVV(N)
             R(K,N)=R(K,N)+A3OVKLN*DDVV(K)
             ROPEN(K,K)=ROPEN(K,K)+A3OVKLN*DOPVV(N)
             ROPEN(K,N)=ROPEN(K,N)+A3OVKLN*DOPVV(K)
         ELSE
C            (KKK)
             R(K,K)=R(K,K)+A3OVKLN*DDVV(K)
             ROPEN(K,K)=ROPEN(K,K)+A3OVKLN*DOPVV(K)
         END IF
 680    CONTINUE
        IF(.NOT. DONE) GO TO 650
C
C
      ELSE
C
C
      CALL THRDIR(IND3OV,HISPIN,.FALSE.,R,ROPEN,ADENS,BDENS,DD,DDVV,
     *            DOPVV,DOP)
C
C
      END IF
C
C       FILL OUT THE TRIANGLE
C
        CALL DSCAL(L3,2.0D+00,R,1)
        CALL DSCAL(L3,2.0D+00,ROPEN,1)
        IF(GOPARR) CALL DDI_GSUMF(1031,R,L3)
        IF(GOPARR) CALL DDI_GSUMF(1032,ROPEN,L3)
C
        CALL CPYSQT(R,S,NUM,1)
        CALL SYMH(S,R,IJINDEX)
        CALL EXPND(S,R,NUM,0)
        IF(NAUXFUN .EQ. 0) CALL DFTTRN2(R,V,S,NUM)
        CALL CPYSQT(ROPEN,S,NUM,1)
        CALL SYMH(S,ROPEN,IJINDEX)
        CALL EXPND(S,ROPEN,NUM,0)
        IF(NAUXFUN .EQ. 0) CALL DFTTRN2(ROPEN,V,S,NUM)
C
C       RESTORE THE DENSITY MATRIX
C
        DO 253 I=1,NUM
           ADENS(I,I)=ADENS(I,I)*2.0D+00
           BDENS(I,I)=BDENS(I,I)*2.0D+00
 253    CONTINUE
      END IF
      DSKWRK = SVDSKW
      END
C*MODULE DFT     *DECK DFTRH4
      SUBROUTINE DFTRH4(ADENS,BDENS,A4OV,R,ROPEN,HISPIN,S,IND4OV,V,
     *                  NAUXFUN,NINTMX)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL HISPIN,DONE,GOPARR,DSKWRK,MASWRK,SVDSKW,DIRSCF,FDIFF
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /IJPAIR/ IJINDEX(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
C
      DIMENSION ADENS(NUM,NUM),BDENS(NUM,NUM),A4OV(*),IND4OV(*),
     *          S(NUM,NUM),R(NUM,NUM),ROPEN(NUM,NUM),V(NUM,NUM)
C
      DONE=.FALSE.
      L3=NUM*NUM
      SVDSKW = DSKWRK
      DSKWRK = .TRUE.
      CALL VCLR(R,1,L3)
      IF(HISPIN) CALL VCLR(ROPEN,1,L3)
C
C     MAKE THE MATRIX REPRESENTATION OF THE DENSITY
C     R(I,J) = <I|(DENSITY)|J>
C
      IF(.NOT. DIRSCF) CALL SEQREW(21)
      IF (.NOT. HISPIN) THEN
C       HALF THE DIAGONAL ELEMENTS
C       SO (IJKK) ISN'T COUNTED TWICE BUT (IJKL) IS
        DO 250 I=1,NUM
           ADENS(I,I)=ADENS(I,I)*0.5D+00
 250    CONTINUE
C
C
        IF(.NOT. DIRSCF) THEN
C
C
 555    CALL PREAD(21,A4OV,IND4OV,ILOOP,NINTMX)
        IF(ILOOP .LE. 0) THEN
           DONE=.TRUE.
           ILOOP = - ILOOP
        END IF
C  F77 MEANS THAT THIS LOOP WILL NOT EXECUTE IF ILOOP=ZERO
        DO 565 I=1,ILOOP
C
C  "IF(N .NE. M) THEN" KEEP US FROM OVER USING THINGS LIKE (2,1,1,1)
C
C   USE LABSIZ TO FIND POSITION (NUM AND BIT SIZE DEPENDANT)
C
                       NPACK = I
                       IF (LABSIZ .EQ. 2) THEN
*I32                     LABEL = IND4OV( 2*NPACK - 1 )
*I32                     IPACK = ISHFT( LABEL, -16 )
*I32                     JPACK = IAND(  LABEL, 65535 )
*I32                     LABEL = IND4OV( 2*NPACK     )
*I32                     KPACK = ISHFT( LABEL, -16 )
*I32                     LPACK = IAND(  LABEL, 65535 )
*I64                     LABEL = IND4OV(NPACK)
*I64                     IPACK = ISHFT( LABEL, -48 )
*I64                     JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                     KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64                     LPACK = IAND( LABEL, 65535 )
                       ELSE IF (LABSIZ .EQ. 1) THEN
*I32                     LABEL = IND4OV(NPACK)
*I32                     IPACK = ISHFT( LABEL, -24 )
*I32                     JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32                     KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I32                     LPACK = IAND( LABEL, 255 )
*I64                     IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                       LABEL = IND4OV( NPACK/2 )
*I64                       IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64                       JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                       LPACK = IAND( LABEL, 255 )
*I64                     ELSE
*I64                       LABEL = IND4OV( (NPACK/2)+1 )
*I64                       IPACK = ISHFT( LABEL, -56 )
*I64                       JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64                       LPACK = IAND( ISHFT( LABEL, -32 ), 255 )
*I64                     END IF
                       END IF
                       K = IPACK
                       L = JPACK
                       N = KPACK
                       M = LPACK
C
             A4OVKLNM=A4OV(I)
C
                 R(K,M)=R(K,M)+A4OVKLNM*ADENS(L,N)
             IF(L .NE. N) THEN
                 R(K,L)=R(K,L)+A4OVKLNM*ADENS(N,M)
                 R(N,M)=R(N,M)+A4OVKLNM*ADENS(K,L)
             END IF
             IF(N .NE. M .AND. L .NE. K) THEN
                 R(L,N)=R(L,N)+A4OVKLNM*ADENS(K,M)
                 R(K,N)=R(K,N)+A4OVKLNM*ADENS(L,M)
                 R(L,M)=R(L,M)+A4OVKLNM*ADENS(K,N)
             ELSE IF(L .NE. K) THEN
                 R(L,M)=R(L,M)+A4OVKLNM*ADENS(K,N)
             ELSE IF(N .NE. M) THEN
                 R(K,N)=R(K,N)+A4OVKLNM*ADENS(L,M)
             END IF
C
C
 565      CONTINUE
      IF( .NOT. DONE) GO TO 555
C
C
      ELSE
C
C
      CALL FORDIR(A4OV,IND4OV,ADENS,BDENS,R,ROPEN,HISPIN)
C
C
      END IF
C
C
C   SYMMETRIZE AND DOUBLE R MATRIX VALUES
C
      CALL DSCAL(L3,2.0D+00,R,1)
      IF(GOPARR) CALL DDI_GSUMF(1033,R,L3)
C
      CALL CPYSQT(R,S,NUM,1)
      CALL SYMH(S,R,IJINDEX)
      CALL EXPND(S,R,NUM,0)
      IF(NAUXFUN .EQ. 0) CALL DFTTRN2(R,V,S,NUM)
C       RESTORE THE DENSITY MATRIX
        DO 251 I=1,NUM
           ADENS(I,I)=ADENS(I,I)*2.0D+00
 251    CONTINUE
      ELSE
C
C
C          OPEN SHELL CASE
C
C
C       HALF THE DIAGONAL ELEMENTS
C       SO (IJKK) ISN'T COUNTED TWICE BUT (IJKL) IS
        DO 252 I=1,NUM
           ADENS(I,I)=ADENS(I,I)*0.5D+00
           BDENS(I,I)=BDENS(I,I)*0.5D+00
 252    CONTINUE
C
C
        IF(.NOT. DIRSCF) THEN
C
C
 655    CALL PREAD(21,A4OV,IND4OV,ILOOP,NINTMX)
        IF(ILOOP .LE. 0) THEN
           DONE=.TRUE.
           ILOOP = - ILOOP
        END IF
C  F77 MEANS THAT THIS LOOP WILL NOT EXECUTE IF ILOOP=ZERO
        DO 665 I=1,ILOOP
C
C  "IF(N .NE. M) THEN" KEEP US FROM OVER USING THINGS LIKE (2,1,1,1)
C
C   USE LABSIZ TO FIND POSITION (NUM AND BIT SIZE DEPENDANT)
C
                       NPACK = I
                       IF (LABSIZ .EQ. 2) THEN
*I32                     LABEL = IND4OV( 2*NPACK - 1 )
*I32                     IPACK = ISHFT( LABEL, -16 )
*I32                     JPACK = IAND(  LABEL, 65535 )
*I32                     LABEL = IND4OV( 2*NPACK     )
*I32                     KPACK = ISHFT( LABEL, -16 )
*I32                     LPACK = IAND(  LABEL, 65535 )
*I64                     LABEL = IND4OV(NPACK)
*I64                     IPACK = ISHFT( LABEL, -48 )
*I64                     JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                     KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64                     LPACK = IAND( LABEL, 65535 )
                       ELSE IF (LABSIZ .EQ. 1) THEN
*I32                     LABEL = IND4OV(NPACK)
*I32                     IPACK = ISHFT( LABEL, -24 )
*I32                     JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32                     KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I32                     LPACK = IAND( LABEL, 255 )
*I64                     IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                       LABEL = IND4OV( NPACK/2 )
*I64                       IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64                       JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                       LPACK = IAND( LABEL, 255 )
*I64                     ELSE
*I64                       LABEL = IND4OV( (NPACK/2)+1 )
*I64                       IPACK = ISHFT( LABEL, -56 )
*I64                       JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64                       LPACK = IAND( ISHFT( LABEL, -32 ), 255 )
*I64                     END IF
                       END IF
                       K = IPACK
                       L = JPACK
                       N = KPACK
                       M = LPACK
C
             A4OVKLNM=A4OV(I)
C
                 R    (K,M)=R    (K,M)+A4OVKLNM*ADENS(L,N)
                 ROPEN(K,M)=ROPEN(K,M)+A4OVKLNM*BDENS(L,N)
             IF(L .NE. N) THEN
                 R    (K,L)=R    (K,L)+A4OVKLNM*ADENS(N,M)
                 ROPEN(K,L)=ROPEN(K,L)+A4OVKLNM*BDENS(N,M)
                 R    (N,M)=R    (N,M)+A4OVKLNM*ADENS(K,L)
                 ROPEN(N,M)=ROPEN(N,M)+A4OVKLNM*BDENS(K,L)
             END IF
             IF(N .NE. M .AND. L .NE. K) THEN
                 R    (L,N)=R    (L,N)+A4OVKLNM*ADENS(K,M)
                 ROPEN(L,N)=ROPEN(L,N)+A4OVKLNM*BDENS(K,M)
                 R    (K,N)=R    (K,N)+A4OVKLNM*ADENS(L,M)
                 ROPEN(K,N)=ROPEN(K,N)+A4OVKLNM*BDENS(L,M)
                 R    (L,M)=R    (L,M)+A4OVKLNM*ADENS(K,N)
                 ROPEN(L,M)=ROPEN(L,M)+A4OVKLNM*BDENS(K,N)
             ELSE IF(L .NE. K) THEN
                 R    (L,M)=R    (L,M)+A4OVKLNM*ADENS(K,N)
                 ROPEN(L,M)=ROPEN(L,M)+A4OVKLNM*BDENS(K,N)
             ELSE IF(N .NE. M) THEN
                 R    (K,N)=R    (K,N)+A4OVKLNM*ADENS(L,M)
                 ROPEN(K,N)=ROPEN(K,N)+A4OVKLNM*BDENS(L,M)
             END IF
C
C
 665      CONTINUE
      IF( .NOT. DONE) GO TO 655
C
      ELSE
         CALL FORDIR(A4OV,IND4OV,ADENS,BDENS,R,ROPEN,HISPIN)
      END IF
C
C   FINISH UP OPEN SHELL CASE...SYMMETRIZE AND DOUBLE R MATRIX VALUES
C
      CALL DSCAL(L3,2.0D+00,R,1)
      CALL DSCAL(L3,2.0D+00,ROPEN,1)
      IF(GOPARR) CALL DDI_GSUMF(1034,R,L3)
      IF(GOPARR) CALL DDI_GSUMF(1035,ROPEN,L3)
C
      CALL CPYSQT(R,S,NUM,1)
      CALL SYMH(S,R,IJINDEX)
      CALL EXPND(S,R,NUM,0)
      IF(NAUXFUN .EQ. 0) CALL DFTTRN2(R,V,S,NUM)
      CALL CPYSQT(ROPEN,S,NUM,1)
      CALL SYMH(S,ROPEN,IJINDEX)
      CALL EXPND(S,ROPEN,NUM,0)
      IF(NAUXFUN .EQ. 0) CALL DFTTRN2(ROPEN,V,S,NUM)
C     RESTORE THE DENSITY MATRIX
      DO 253 I=1,NUM
         ADENS(I,I)=ADENS(I,I)*2.0D+00
         BDENS(I,I)=BDENS(I,I)*2.0D+00
 253  CONTINUE
      END IF
C          END OF OPEN SHELL CASE
C
      DSKWRK = SVDSKW
      RETURN
      END
C*MODULE DFT     *DECK DFTTRN2
      SUBROUTINE DFTTRN2(A,B,C,N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(N,N),B(N,N),C(N,N)
C
C     THIS ROUTINE TRANSFORMS A BY A=(B DAGGER)*A*B
C     CONSIDER REPLACING WITH TFSQU TO GET PARALLEL, BUT
C     LOSE BLAS?
C
      CALL MRARBR(A,N,N,N,B,N,N,C,N)
      CALL MRTRBR(B,N,N,N,C,N,N,A,N)
      END
C*MODULE DFT     *DECK DFTWRIT
      SUBROUTINE DFTWRIT(R,ROPEN,RSA,RSB,D,VS,V,
     *                   AFOCK,BFOCK,S,HISPIN,ADENS,BDENS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL HISPIN,DIRSCF,FDIFF
C
      PARAMETER (MXATM=2000)
C
      COMMON /CONV  / DENTOL,EN,ETOT,EHF,EHF0,DIFF,ITER,ICALP,ICBET
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OPTSCF/ DIRSCF,FDIFF
C
      DIMENSION R    (NUM+NAUXFUN,NUM+NAUXFUN),
     *          ROPEN(NUM+NAUXFUN,NUM+NAUXFUN),
     *          RSA  (NUM+NAUXFUN,NUM+NAUXFUN),
     *          RSB  (NUM+NAUXFUN,NUM+NAUXFUN),
     *          D    (NUM+NAUXFUN,NUM+NAUXFUN),
     *          VS   (NUM+NAUXFUN,NUM+NAUXFUN),
     *          V    (NUM+NAUXFUN,NUM+NAUXFUN),
     *          AFOCK(NUM*(NUM+1)/2),
     *          BFOCK(NUM*(NUM+1)/2),
     *          S(NUM+NAUXFUN,NUM+NAUXFUN),
     *          ADENS(NUM,NUM),BDENS(NUM,NUM)
C
      L2 = (NUM*NUM+NUM)/2
      L4 = NUM+NAUXFUN
      L5=(L4+L4*L4)/2
      L6 = L4*L4
C
C     TAKE INTO ACCOUNT THE DIFFERENCE OF THE FOCK FROM THE ENERGY
C     HARTREE FOCK CODE WILL INCORRECTLY DO R*DENSITY LATER ON
C     THIS SUBTRACTS THAT BACKOUT
C
C     COMPUTE DFT CONTRIBUTION TO THE TOTAL ENERGY
C
C     IN THE MO BASIS DENSITY MATRIX IS VERY SIMPLE (DIAGONAL)
C
      EXENA=0.0D+00
      EXENB=0.0D+00
      EXENC=0.0D+00
C
C     THE OVERLAP MATRIX--MAKE VS
C
      IF(NAUXFUN .EQ. 0) THEN
        CALL DAREAD(IDAF,IODA,S,L2,12,0)
      ELSE
        CALL DAREAD(IDAF,IODA,S,L5,350,0)
      END IF
      CALL CPYTSQ(S,D,L4,1)
      CALL MRARBR(D,L4,L4,L4,V,L4,L4,VS,L4)
C
C    TRANSFORM TO AO BASIS BY RAO=(VS)RMO(VS DAGGER)
C
      CALL MRARTR(R,L4,L4,L4,VS,L4,L4,S,L4)
      CALL MRARBR(VS,L4,L4,L4,S,L4,L4,R,L4)
      IF(HISPIN) THEN
        CALL MRARTR(ROPEN,L4,L4,L4,VS,L4,L4,S,L4)
        CALL MRARBR(VS,L4,L4,L4,S,L4,L4,ROPEN,L4)
      END IF
C
      IF(HISPIN) THEN
          CALL MRARTR(RSA,L4,L4,L4,VS,L4,L4,S,L4)
          CALL MRARBR(VS,L4,L4,L4,S,L4,L4,RSA,L4)
          CALL MRARTR(RSB,L4,L4,L4,VS,L4,L4,S,L4)
          CALL MRARBR(VS,L4,L4,L4,S,L4,L4,RSB,L4)
          DO I=1,NUM
            DO J=1,NUM
              EXENA=EXENA+RSA  (I,J)*ADENS(I,J)
              EXENC=EXENC+R    (I,J)*ADENS(I,J)
              EXENC=EXENC+ROPEN(I,J)*BDENS(I,J)
              EXENB=EXENB+RSB  (I,J)*BDENS(I,J)
            END DO
          END DO
          EXENC=-EXENC*0.5D+00
      ELSE
          CALL MRARTR(RSA,L4,L4,L4,VS,L4,L4,S,L4)
          CALL MRARBR(VS,L4,L4,L4,S,L4,L4,RSA,L4)
          DO I=1,NUM
            DO J=1,NUM
              EXENA=EXENA+RSA(I,J)*ADENS(I,J)
              EXENC=EXENC+R  (I,J)*ADENS(I,J)
            END DO
          END DO
          EXENB=EXENA
          EXENC=-EXENC
      END IF
C
C     ADD DFT MATRIX R (AND ROPEN) TO FOCK MATRIX
C
      IF (.NOT. HISPIN) THEN
C
C     HANDLE DIRECT RUNS, BY SUBTRACTING OUT PREVIOUS FOCK CONTRIBUTION
C
        IF(DIRSCF .AND. FDIFF .AND. ITER.NE.1) THEN
         CALL DAREAD(IDAF,IODA,D,L6,340,0)
         CALL CPYSQT(D,S,L4,1)
         CALL VSUB(S,1,AFOCK,1,AFOCK,1,L2)
        END IF
        CALL CPYSQT(R,S,L4,1)
        CALL VADD (AFOCK,1,S,1,AFOCK,1,L2)
      ELSE
C
C     HANDLE DIRECT RUNS, BY SUBTRACTING OUT PREVIOUS FOCK CONTRIBUTION
C
        IF(DIRSCF .AND. FDIFF .AND. ITER.NE.1) THEN
         CALL DAREAD(IDAF,IODA,D,L6,340,0)
         CALL CPYSQT(D,S,L4,1)
         CALL VSUB(S,1,AFOCK,1,AFOCK,1,L2)
         CALL DAREAD(IDAF,IODA,D,L6,341,0)
         CALL CPYSQT(D,S,L4,1)
         CALL VSUB(S,1,BFOCK,1,BFOCK,1,L2)
        END IF
        CALL CPYSQT(R,S,L4,1)
        CALL VADD(AFOCK,1,S,1,AFOCK,1,L2)
        CALL CPYSQT(ROPEN,S,L4,1)
        CALL VADD(BFOCK,1,S,1,BFOCK,1,L2)
      END IF
C
C     SAVE FOCK MATRIX FOR GRADIENTS, DIRECT RUNS ETC.
C
      CALL DAWRIT(IDAF,IODA,R,L6,340,0)
      IF(HISPIN) CALL DAWRIT(IDAF,IODA,ROPEN,L6,341,0)
      END
