C 20 AUG 07 - MC,DGF - TDDFT GRADIENT CHANGES, SYNCH FMO COMMONS
C 20 AUG 07 - DDK,LVS - VNNDER: CALLING ARGS TO GEXREP, EWALD SUMS
C 24 MAR 07 - HM  - ADD CALLS FOR MCP ANALYTIC GRADIENT (SPDFG AO'S)
C 14 FEB 07 - MWS - VNNDER: PROVIDE DYNAMIC STORAGE FOR EFP TEMPS
C 22 DEC 06 - DGF,TN - STANDARDISE FMO LAGRANGIAN, FMO/EFP GRADIENT
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 29 MAR 06 - FLG,JK - STVDER: LAGRANGIAN FOR ELONGATION METHOD
C 13 MAR 06 - MWS - SEPARATE DAF RECORD NOS OF RHF/ROHF MP2 FROM UHF MP2
C                   SEPARATE LAGRANGIAN DAF RECORD NOS FOR GVB,MCSCF,CI
C 17 JAN 06 - HL  - VNNDER: MODIFY GRADIENTS FOR MD
C 14 NOV 05 - DGF - SMALL TOUCH FOR FMO/PCM
C 19 SEP 05 - MWS - VNNDER: PERFORM THE EFP GLOBAL SUMS ONLY IF EFP
C 13 JUL 05 - MWS - VNNDER: ADDITIONAL DYN.STORE FOR GCTANA
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  1 JUN 05 - CMA - CHANGES FOR ZAPT2 GRADIENT
C  1 JUN 05 - MWS - VNNDER: IMPROVE EFP DEBUGGING PRINT
C 30 APR 05 - DGF - CHANGES FOR FMO-MCSCF GRADIENT
C 14 MAR 05 - HL  - VNNDER: EFP CHARGE TRANSFER, EX-REP GRADIENTS
C  7 MAR 05 - IA  - VNNDER: READ DYN.POL. CONTRIBUTION FROM DAF
C 22 FEB 05 - MWS - VNNDER: JUST BEAUTIFICATION OF DEBUG LINES
C 13 FEB 05 - HMN - VNNDER: EFP PARALLELIZATION AND MD HOOKS
C 13 FEB 05 - JS  - VNNDER: EFP MP2 CHANGES
C  5 FEB 05 - MWS - PAD COMMON FMORUN AND COMMON MCINP
C 23 JUL 04 - DGF - EIJDEN: READ CORRECT CI LAGRANGIAN IF SPHER.HARM.
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C  9 DEC 03 - OQ  - SYNCH ARGUMENTS TO DIPDER
C  4 NOV 03 - TN  - STVDER: CHANGES FOR DOUGLAS-KROLL
C  3 SEP 03 - SPW - EIJDEN: CHANGES FOR CIS GRADIENT
C  3 JUL 03 - JMM - SUPPRESS PRINTING FOR MONTE CARLO JOBS
C 14 JAN 03 - JI  - ADD HOOK FOR ORMAS CI STEP
C  7 AUG 02 - HT  - USE DYNAMIC MEMORY FOR EFP PAULI REPULSION
C 17 APR 02 - MWS - SYNCH UP FRGINF COMMON
C 24 JAN 02 - CMA - EIJDEN: ALLOW FOR UHF-MP2
C 21 DEC 01 - MWS - DENDD1: RETURN TRUE SCF DENSITY FOR SCF STAGE OF CI
C  1 AUG 01 - JI  - TWEAKS FOR GENERAL CI
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 13 JUN 00 - DGF,TN - IMPLEMENT INTERNALLY UNCONTRACTED RESC METHOD
C 20 FEB 01 - MWS - PAD PAULMO COMMON
C 11 OCT 00 - PB  - VNNDER: INTERFACING EFP AND PCM
C 10 APR 00 - MWS - REMOVE STATIC MEMORY FROM COSMO COMMONS
C 25 MAR 00 - DGF - CODE NESC GRADIENT
C  8 MAR 00 - KKB/LNB - COSMO CONTRIBUTION TO GRADIENT
C 16 FEB 00 - VK  - STVDER,CIGRDM: CHANGES FOR PRESENCE OF BUFFER MO-S
C 14 JAN 00 - VK  - EIJDEN: SPECIAL TREATMENT OF RHF MO FREEZING
C 21 DEC 99 - MWS - MAKE GUGWFN COMMON CONSISTENT
C 29 FEB 00 - DGF,TN - CODE RESC GRADIENT
C 21 DEC 99 - MWS - MAKE GUGWFN COMMON CONSISTENT
C  1 DEC 98 - BMB - UPDATE ECPGRD CALLS FOR NEW ECP CODE
C 12 NOV 98 - GDF - SYMEG: DELETE USE OF ISOOUT
C 26 OCT 98 - MAF - EIJDEN: ALLOW FOR USE OF SPHERICAL HARMONICS
C 12 APR 98 - MWS - EIJDEN: CHANGES FOR DETERMINANT DENSITY MATRIX
C 27 FEB 98 - MWS - DENDD1: ALLOW FOR CI OTHER THAN GUGA
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C 16 JUL 97 - GNM - VNNDER: SKIP SOME TERMS FOR TIP*P POTENTIALS
C 14 FEB 97 - MWS - STVDER: PRINT LAGRANGIAN AND DENSITY IF DEBUGGING
C 18 DEC 96 - JHJ - STVDER: CALLS TO DPAULA, QMADD, DIPDER CALL CHANGED.
C                   VNNDER: CALL TO DPAULR.
C 26 NOV 96 - SPW - EIJDEN: INCLUDE FROZEN CORES IN CI LAGRANGIAN TRANSF
C 17 OCT 96 - SPW - EIJDEN,STVDER: CHANGES FOR CI GRADIENTS, ADD CIGRDM
C 18 SEP 96 - MWS - VNNDER: PASS DYNAMIC MEMORY TO DREPUL
C 13 JUN 96 - VAG - CHANGES TO INTRODUCE CITYP VARIABLE
C  5 MAR 96 - MWS - CHANGE SHELL SYMMETRY PACKING COMMONS
C 11 NOV 95 - SPW - EIJDEN: DAREAD ENERGY WEIGHTED MP2 DENSITY MATRIX
C 27 MAR 95 - SPW - CHANGE ARGUMENTS TO DIPDER CALL
C 29 DEC 94 - TLW - STVDER,VNNDER: INCLUDE ELECTRIC FIELD TERMS
C 17 NOV 94 - WC  - REDIMENSION THE FGRAD COMMON
C 12 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 27 MAY 94 - PND - STVDER,VNNDER: INCLUDE FRAGMENT-FRAGMENT GRADIENT
C  4 APR 94 - MWS - STVDER: SEPARATE ZRF AND EFP CALLS
C  4 NOV 93 - MH  - TVDER: CHARMM EXTERNAL CHARGE PERT. ADDED
C 10 AUG 93 - MWS - SYMEG BROUGHT TO THIS FILE
C  5 AUG 93 - BMB - ADDED F AND G GRADIENT, CHANGED DYNAMIC MEMORY
C                   CALLS SO STVDER HANDLES ALL DYNAMIC MEMORY AND I/O
C  1 APR 93 - PND - ADDITIONAL EF INTERFACE CHANGES
C  4 MAR 93 - JHJ - COMPUTE SCRF GRADIENT CONTRIBUTION
C 18 DEC 92 - MWS - SDER: REMOVE FAC COMPUTATION FROM J SHELL LOOP
C  5 NOV 92 - TLW - MOVE NEXT VALUE BALANCING
C 14 MAY 92 - MWS - FIX NEXT VALUE BALANCING IN SDER
C  5 MAY 92 - PND - SGRAD, FRAGMENT GRADIENT CONTRIBUTIONS
C  5 MAY 92 - PND - STVDER, FRAGMENT GRADIENT CONTRIBUTIONS
C 23 MAR 92 - MWS - MOVE TVDER BEFORE ECP TERMS, PRINT ECP TIME
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C  4 JAN 92 - TLW - MAKE WRITES PARALLEL;ADD COMMON PAR
C 24 OCT 91 - JHJ - INTRODUCED ICALC IN COMMON /ZRFPAR/.
C 23 JUL 91 - JHJ - TVDER:MOVE ZRFDER AND EFXDER CALLS TO FRGDER.
C                   STVDER:CALL FRGDER. MOVED NUC. EFC GRAD. TO DNUCC
C                   AND EFC GRADIENTS TO EFCDER.REMOVED CHECK FROM
C                   ECPDER CALL SINCE THAT IS DONE IN ECPDER.
C 26 FEB 91 - DRG - TVDER:CALL ZRFDER (FIXED FIELD DERIVATIVES).
C  1 NOV 90 - MWS - GLIMIT: ALLOW FOR GHOST ATOMS.
C 12 SEP 90 - MWS - CHECK RUNS AVOID ECP GRADIENT INTEGRALS
C 15 AUG 90 - TLW - ADD H7(7) TO COMMON HERMIT AND W7(7) TO COMMON
C                   WERMIT.  ADD 7TH ELEMENT TO MIN AND MAX ARRAYS.
C 15 AUG 90 - DRG - FIXED INDUCED MOMENT DERIVATIVES FOR OPTIMIZATIONS.
C                   CALL EFPDER IN TVDER AND DNUCP IN SGRAD.
C 24 MAY 90 - DRG - OCTUPOLE DERIVATIVES:CALL EFODER IN TVDER AND
C                   DNUCO IN SGRAD.
C 23 MAY 90 - DRG - QUADRUPOLE DERIVATIVES:CALL EFQDER IN TVDER AND
C                   DNUCQ IN SGRAD.
C 17 MAY 90 - DRG - DIPOLE DERIVATIVES:CALL EFDDER IN TVDER AND DNUCD
C                   IN SGRAD.
C 14 MAY 90 - DRG - NUCLEAR GRADIENTS WRT EFC'S ADDED TO SGRAD.
C                   EFECTIVE POINT CHARGE GRADIENTS ADDED IN TVDER.
C  7 AUG 90 - TLW - ADD CF AND CG TO COMMON NSHEL
C 29 NOV 89 - MWS - ADD NECP TO /INTOPT/
C 16 SEP 89 - KAN - STVDER: ECP RUNS CALL ECPDER
C  6 MAY 89 - MWS - USE PRTRI TO PRINT LAGRANGIAN AFTER EIJDEN
C 20 MAR 89 - MWS - SIMPLIFY A BIT IN EIJDEN AND SGRAD
C  2 MAR 89 - MWS - FIX CALL TO EGOUT IN SGRAD
C 25 FEB 89 - STE - SGRAD: DEFINE EXETYP
C 16 JAN 89 - MWS - SUPPRESS PRINT OF HELLMANN-FEYNMAN FORCE,
C                   RENAME MODULE FROM INTG1 TO GRD1,
C                   USE DAWRIT RATHER THAN WRTGRD.
C 20 DEC 88 - MWS - USE EXETYP=INTG1 FOR MODEST AMOUNTS OF DEBUG
C 10 AUG 88 - MWS - MXSH,MXGSH,MXGTOT FROM 120,10,440 TO 1000,30,5000
C 24 MAY 88 - MWS - USE PARAMETERS TO DIMENSION COMMONS
C  8 MAY 88 - MWS - UP AO-S FROM 255 TO 2047
C 22 APR 88 - MWS - USE /GVBWFN/ IN EIJDEN
C 14 APR 88 - MWS - EIJDEN CHANGED FOR NEW ROHF CODE
C 14 FEB 88 - MWS - CHANGE RYS ROOTS TO 9 ROOTS
C 14 NOV 87 - STE - USE EXETYP
C  8 MAY 87 - STE - USE PARAMETERS; VINT: PROTECT AGAINST UNDERFLOWS
C  7 JUL 86 - MWS - SANITIZE FLOATING POINT CONSTANTS
C 31 OCT 85 - STE - USE GENERIC EXP,SQRT; EIJDEN: USE DCOPY,VCLR
C                   TVDER: PRESERVE RESTART VALUE
C  2 APR 84 - STE - REPLACE DV123 WITH RT123
C  2 OCT 82 - MWS - CONVERT FOR THE IBM MACHINE
C
C*MODULE GRD1    *DECK CIGRDM
      SUBROUTINE CIGRDM
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL SOME,GOPARR,DSKWRK,MASWRK,EFLDL
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXAO=8192,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /EFLDC / EVEC(3),EFLDL
      COMMON /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      COMMON /EFPOTD/ METHOF,ISCCHG,MOVE
      COMMON /FMCOM / X(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
      DATA ZRF/8HZRF     /
C
      IPOT=IEFC+IEFD+IEFQ+IEFO+IEFP+IREP
C
      L1 = NUM
      L2 = (NUM*NUM+NUM)/2
      L3 = NUM*NUM
C
      SOME = NPRINT.NE.-5 .AND. MASWRK
      T0 = 0.0D+00
      CALL TSECND(T0)
C
C           STORAGE FOR MCP CAN BE TRUNCATED TO TRIVIAL AMOUNTS.
C           THE 900'S ARE FIRMLY EMBEDDED IN THE MCP CODE, SO YOU'D
C           NEED TO DO MORE THAN JUST CHANGE THE TWO SIZES HERE.
      IF(IECP.EQ.5) THEN
         NATMCP=NAT
         MXVO  =900
         MXCO  =900
      ELSE
         NATMCP=1
         MXVO=1
         MXCO=1
      END IF
C
      CALL VALFM(LOADFM)
      LDA    = LOADFM + 1
      LDB    = LDA + L3
      LVEC   = LDB + L3
      LDBF   = LVEC + L3
      LWRK   = LDBF + L2
      LOCC   = LWRK + L2
      LDD    = LOCC + L1
      LOV    = LDD  + 3*NATMCP
      LDV    = LOV  + MXCO*MXVO
      LAST   = LDV  + MXCO*MXVO*3
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      CALL DENDD1(X(LDA),X(LDB),L2)
C
C     ----- EFFECTIVE CORE POTENTIAL MODIFICATIONS -----
C     DYNAMIC MEMORY FOR ECDINT FOLLOWS PREVIOUS STORAGE
C
C     *** MODEL CORE POTENTIAL POTENTIAL MODIFICATIONS WERE
C     *** ADDED BY HIROTOSHI MORI & EISAKU MIYOSHI
C
      IF(IECP.GT.0) THEN
         IF (IECP.EQ.5) THEN
            CALL MCPGRD(DE, X(LDA), X(LDD), X(LOV), X(LDV), MXCO, MXVO)
         ELSE
            CALL ECPGRD(DE, X(LDA), EXETYP, SOME)
         END IF
         CALL TSECND(T0)
      END IF
C
C     ----- GET SELF-CONSISTENT REACTION FIELD GRADIENT -----
C     THIS SOLVENT OPTION ISN'T ENABLED FOR GRADIENTS YET
C
      IF(IZRF.EQ.1) THEN
         CALL ABRT
         DUMMY = 0.0D+00
         CALL DIPDER(X(LDA),DUMMY,DUMMY,DUMMY,L2,DUMMY,0,0,0,ZRF)
         IF(SOME) THEN
            CALL TSECND(T1)
            TG = T1-T0
            WRITE(IW,9040) TG
            T0 = T1
         END IF
      END IF
C
C     ----- EFFECTIVE FRAGMENT GRADIENT -----
C     EFFT CALCULATES TOTAL FORCE AND TORQUE ON EACH FRAGMENT
C
      IF(IPOT.GT.0) THEN
         CALL ABRT
C    BUFFER MO DENSITY SHOULD BE SUBTRACTED FROM THE TOTAL
         IF(NBUFMO.GT.0)THEN
           CALL DCOPY(L2,X(LDA),1,X(LWRK),1)
           CALL DENDBF(X(LDA),X(LDB),X(LDBF),X(LVEC),X(LOCC),NBUFMO,
     *                 L1,L2,L3)
         END IF
         CALL FRGDER(X(LDA),L2)
         IF(NBUFMO.GT.0) CALL DCOPY(L2,X(LWRK),1,X(LDA),1)
         IF(MOVE.NE.0) CALL EFFT
         IF(SOME) THEN
            CALL TSECND(T1)
            TG = T1-T0
            WRITE(IW,9050) TG
            T0 = T1
         END IF
      END IF
C
C     ----- GET ELECTRIC FIELD CONTRIBUTION -----
C
      IF (EFLDL) THEN
         CALL EFLD2(X(LDA),L2)
         IF(SOME) THEN
            CALL TSECND(T1)
            TG = T1-T0
            WRITE(IW,9060) TG
            T0 = T1
         END IF
      END IF
C
      CALL RETFM(NEED)
C
C     ----- GATHER ALL GRADIENT CONTRIBUTIONS -----
C     ----- SAVE COMPLETED 1E-GRADIENT -----
C
      IF(GOPARR) CALL DDI_GSUMF(1500,DE,3*NAT)
      CALL DAWRIT(IDAF,IODA,DE,3*NAT,3,0)
C
C     ----- CHECK TIME -----
C
      IREST=4
      IST=1
      JST=1
      KST=1
      LST=1
      CALL TEXIT(1,IREST)
      RETURN
C
 9040 FORMAT(1X,'TIME TO DO     SCRF GRADIENT INTEGRALS=',F10.2)
 9050 FORMAT(1X,'TIME TO DO FRAGMENT GRADIENT INTEGRALS=',F10.2)
 9060 FORMAT(1X,'TIME TO DO ELECTRIC FIELD CONTRIBUTION=',F10.2)
      END
C*MODULE GRD1    *DECK DENDD1
      SUBROUTINE DENDD1(DA,DB,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PSILVL/ IPSI,ISKPRP
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DIMENSION DA(L2),DB(L2)
C
      DATA RHF,RMC/8HRHF     ,8HMCSCF   /
      DATA RNONE/8HNONE    /
C
C        READ TOTAL DENSITY MATRIX
C        -DENDD1- DIFFERS FROM -WFNDEN- IN THAT IT ADDS DA+DB.
C
      CALL DAREAD(IDAF,IODA,DA,L2,16,0)
      IF(SCFTYP.EQ.RHF) RETURN
      IF(SCFTYP.EQ.RMC) RETURN
      IF(CITYP.NE.RNONE  .AND.  IPSI.EQ.1) RETURN
C
      CALL DAREAD(IDAF,IODA,DB,L2,20,0)
C
C        RETURN TOTAL DENSITY IN -DA-
C
      CALL VADD(DA,1,DB,1,DA,1,L2)
      RETURN
      END
C*MODULE GRD1    *DECK DERI
      SUBROUTINE DERI(DXDI,DYDI,DZDI,X,Y,Z,LIT,LJT,AI)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION  X  (6,1), Y  (6,1), Z  (6,1)
      DIMENSION DXDI(5,1),DYDI(5,1),DZDI(5,1)
C
C     ----- DXDI ... -----
C
      DO 110 J=1,LJT
        DXDI(1,J)=X(2,J)*(AI+AI)
        DYDI(1,J)=Y(2,J)*(AI+AI)
        DZDI(1,J)=Z(2,J)*(AI+AI)
  110 CONTINUE
C
      IF(LIT.EQ.1) RETURN
C
      DO 120 I=2,LIT
        DO 120 J=1,LJT
          DXDI(I,J)=X(I+1,J)*(AI+AI) - X(I-1,J)*(I-1)
          DYDI(I,J)=Y(I+1,J)*(AI+AI) - Y(I-1,J)*(I-1)
          DZDI(I,J)=Z(I+1,J)*(AI+AI) - Z(I-1,J)*(I-1)
  120 CONTINUE
C
      RETURN
      END
C*MODULE GRD1    *DECK DTXYZ
      SUBROUTINE DTXYZ(XT,YT,ZT,XS,YS,ZS,NI,NJ,AJ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION XT(6,*),YT(6,*),ZT(6,*)
      DIMENSION XS(6,*),YS(6,*),ZS(6,*)
      PARAMETER (THREE=3.0D+00)
C
      DO 10 I=1,NI
        XT(I,1)=(XS(I,1  )       - XS(I,3  )*(AJ+AJ))*AJ
        YT(I,1)=(YS(I,1  )       - YS(I,3  )*(AJ+AJ))*AJ
        ZT(I,1)=(ZS(I,1  )       - ZS(I,3  )*(AJ+AJ))*AJ
   10 CONTINUE
C
      IF(NJ.EQ.1) RETURN
C
      DO 20 I=1,NI
        XT(I,2)=(XS(I,2  )*THREE - XS(I,4  )*(AJ+AJ))*AJ
        YT(I,2)=(YS(I,2  )*THREE - YS(I,4  )*(AJ+AJ))*AJ
        ZT(I,2)=(ZS(I,2  )*THREE - ZS(I,4  )*(AJ+AJ))*AJ
   20 CONTINUE
C
      IF(NJ.EQ.2) RETURN
C
      DO 30 J=3,NJ
        FACT1 = J+J-1
        FACT2 = (J-1)*(J-2)
        FACT2 = FACT2/2
        DO 30 I=1,NI
          XT(I,J)=(XS(I,J  )*FACT1 - XS(I,J+2)*(AJ+AJ))*AJ
     *                             - XS(I,J-2)*FACT2
          YT(I,J)=(YS(I,J  )*FACT1 - YS(I,J+2)*(AJ+AJ))*AJ
     *                             - YS(I,J-2)*FACT2
          ZT(I,J)=(ZS(I,J  )*FACT1 - ZS(I,J+2)*(AJ+AJ))*AJ
     *                             - ZS(I,J-2)*FACT2
   30 CONTINUE
C
      RETURN
      END
C*MODULE GRD1    *DECK DVINT
      SUBROUTINE DVINT
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /DSTV  / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,
     *                XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      COMMON /HERMIT/ H1,H2(2),H3(3),H4(4),H5(5),H6(6),H7(7),
     *                H8(8),H9(9),H10(10)
      COMMON /WERMIT/ W1,W2(2),W3(3),W4(4),W5(5),W6(6),W7(7),
     *                W8(8),W9(9),W10(10)
C
      DIMENSION H(28),W(28),MIN(7),MAX(7)
C
      EQUIVALENCE (H(1),H1),(W(1),W1)
C
      DATA MIN /1,2,4,7,11,16,22/
      DATA MAX /1,3,6,10,15,21,28/
C
C     ----- GAUSS-HERMITE QUADRATURE USING MINIMUM POINT FORMULA -----
C
      XINT = ZERO
      YINT = ZERO
      ZINT = ZERO
      NPTS = (NI+NJ+1-2)/2+1
      IMIN = MIN(NPTS)
      IMAX = MAX(NPTS)
      DO 340 I = IMIN,IMAX
         DUM = H(I)*T
         PTX = DUM+X0
         PTY = DUM+Y0
         PTZ = DUM+Z0
         PX = PTX-CX
         PY = PTY-CY
         PZ = PTZ-CZ
         AX = PTX-XI
         AY = PTY-YI
         AZ = PTZ-ZI
         BX = PTX-XJ
         BY = PTY-YJ
         BZ = PTZ-ZJ
         GO TO (180,170,160,140,120,100),NI
C
  100    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
  120    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
  140    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
  160    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
  170    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
  180    GO TO (320,300,280,260,240,220,200),NJ
C
  200    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  220    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  240    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  260    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  280    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  300    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  320    DUM = W(I)
         XINT = XINT+DUM*PX
         YINT = YINT+DUM*PY
         ZINT = ZINT+DUM*PZ
  340 CONTINUE
      RETURN
      END
C*MODULE GRD1    *DECK EIJDEN
      SUBROUTINE EIJDEN(EPS,V,E,IA,WRK,L1,L2,L3,L0,ECI)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION METHOD
C
      LOGICAL CANONC,FCORE,FORS,NOCI,EKT,LINSER,MFRZ
C
      PARAMETER (MXATM=2000, MXAO=8192, MXRT=100, MXNORO=250)
C
      DIMENSION V(L1,L0),EPS(L2),E(*),IA(L1),WRK(L1),ECI(L0,L0)
C
C        E(L1) FOR RHF, E(L3) FOR GVB/MCSCF/CI
C        NOTE FOR CI RUNS, -ECI- AND -E- MUST BE EQUIVALENT BY CALL
C
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,SDET,SZDET,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORBDT,
     *                NADET,NBDET,KDET,KSTDET,IROOT,IPURES,MAXW1,NITDET,
     *                MAXP,NCIDET,IGPDET,KSTSYM,NFTGCI
      COMMON /FMOOPT/ ESPSCA(9),RESPAP(2),RESPPC(2),RESDIM,RESTRI(4),
     *                RCORSD,RESPCT,CONVFG,CNVDMP,COROFF,IXESP,MXITFG,
     *                NGUESS,NBSSE,MODORB,MODPAR,IRSTSTP,IRSTLAY,NPRFMO,
     *                NFMOPAL,MODPRP,MAXL1C,IPIEDA,MODGRD,MODESP,IVMUL
      COMMON /GUGWFN/ NFZC,NMCC,NDOC,NAOS,NBOS,NALP,NVAL,NEXT,NFZV,
     *                IFORS,IEXCIT,ICICI,NOIRR
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NMOGVB,NCONF(MXAO),NHAM
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,
     *                NORBMC,NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /MFRPAR/ MFRZ,NUMFRZ,IFRZ(MXAO)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
C
      PARAMETER (ZERO=0.0D+00, HALF=0.5D+00)
C
      DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF    /
      DATA  GVB,RMC/8HGVB     ,8HMCSCF   /
      DATA ALDET,GUGA,GENCI,ORMAS
     *      /8HALDET   ,8HGUGA    ,8HGENCI   ,8HORMAS   /
      DATA CIS/8HCIS     /
C
C     ----- COMPUTE "ENERGY WEIGHTED DENSITY MATRIX" -----
C     THIS QUANTITY IS ACTUALLY THE LAGRANGIAN MATRIX,
C     BACKTRANSFORMED INTO THE AO BASIS.
C
C           RHF CASE
C
      IF(SCFTYP.EQ.RHF .AND. MPLEVL.EQ.0 .AND.
     *   CITYP.NE.GUGA .AND. CITYP.NE.CIS) THEN
         IF(MFRZ) THEN
           L2 = (L1*L1+L1)/2
           CALL DAREAD(IDAF,IODA,EPS,L2,36,0)
           RETURN
         END IF
         CALL DAREAD(IDAF,IODA,V,L3,15,0)
         CALL DAREAD(IDAF,IODA,E,L1,17,0)
         IJ = 0
         DO 130 I = 1,L1
            DO 120 J = 1,I
               IJ = IJ+1
               DUM = ZERO
               DO 100 K = 1,NA
                  DUM = DUM-E(K)*V(I,K)*V(J,K)
  100          CONTINUE
               EPS(IJ) = DUM+DUM
  120       CONTINUE
  130    CONTINUE
         IF(NFG.NE.0) THEN
           L2 = (L1*L1+L1)/2
           IF(IAND(MODGRD,1).EQ.0) CALL FMOLAG(EPS,V,WRK,L1,L2)
C          V IS USED AS TEMPORARY STORAGE
         ENDIF
         RETURN
      END IF
C
C           UHF AND ROHF CASE
C
      IF((SCFTYP.EQ.UHF .OR. SCFTYP.EQ.ROHF)  .AND.  MPLEVL.EQ.0) THEN
         L2 = (L1*L1+L1)/2
         CALL DAREAD(IDAF,IODA,EPS,L2,36,0)
         IJ = 0
         DO 160 I = 1,L1
            IJ = IJ+I
            EPS(IJ) = HALF*EPS(IJ)
  160    CONTINUE
         RETURN
      END IF
C
C     CLOSED SHELL CIS CASE: JUST READ IT FROM DAF
C
      IF(SCFTYP.EQ.RHF .AND. CITYP.EQ.CIS) THEN
         L2 = (L1*L1+L1)/2
         CALL DAREAD(IDAF,IODA,EPS,L2,309,0)
         RETURN
      END IF
C
C     CLOSED SHELL MP2 CASE: JUST READ IT FROM DAF
C
      IF(SCFTYP.EQ.RHF .AND. MPLEVL.EQ.2) THEN
         L2 = (L1*L1+L1)/2
         CALL DAREAD(IDAF,IODA,EPS,L2,309,0)
C        THE OLD MP2 WAS NOT TO ADJUST LAGRANGIAN.
         IF(NFG.NE.0.AND.IAND(MODGRD,1).EQ.0)
     *     CALL FMOLAG(EPS,V,WRK,L1,L2)
         RETURN
      END IF
C
C     OPEN SHELL ZAPT2 CASE: JUST READ IT FROM DAF
C     OPEN SHELL OSPT=RMP2 SHOULD NOT GET THIS FAR (INPUTA ABRT)
C
      IF(SCFTYP.EQ.ROHF .AND. MPLEVL.EQ.2) THEN
         L2 = (L1*L1+L1)/2
         CALL DAREAD(IDAF,IODA,EPS,L2,309,0)
         RETURN
      END IF
C
C     UMP2 CASE: READ ALPHA AND BETA ENERGY WEIGHTED DENSITY, AND ADD
C
      IF(SCFTYP.EQ.UHF .AND. MPLEVL.EQ.2) THEN
         L2 = (L1*L1+L1)/2
         CALL DAREAD(IDAF,IODA,EPS,L2,419,0)
         CALL DAREAD(IDAF,IODA,E  ,L2,429,0)
         DO N=1,L2
            EPS(N)=EPS(N)+E(N)
         ENDDO
         RETURN
      END IF
C
C           ONLY GVB/MCSCF/CI GET THIS FAR...
C     CIS DOES NOT PASS THROUGH HERE, OTHERWISE ONLY GUGA HAS GRADIENTS.
C
      DO I = 1,L1
         IA(I) = (I*I-I)/2
      ENDDO
C
      IF(SCFTYP.EQ.GVB) NORB = NMOGVB
      IF(SCFTYP.EQ.RMC) THEN
         IF(CISTEP.EQ.ALDET) NORB = NCORSV+NACT
         IF(CISTEP.EQ.GENCI) NORB = NCORSV+NACT
         IF(CISTEP.EQ.ORMAS) NORB = NCORSV+NACT
         IF(CISTEP.EQ.GUGA)  NORB = NMCC+NDOC+NAOS+NBOS+NALP+NVAL
      END IF
      IF(CITYP.EQ.GUGA) NORB = NFZC+NDOC+NAOS+NBOS+NALP+NVAL
C
C     ----- READ MOS AND LAGRANGIAN FROM DISK -----
C
      CALL DAREAD(IDAF,IODA,V,L3,15,0)
C
                        L3ORB = NORB*NORB
      IF(SCFTYP.EQ.RMC) L3ORB = L3
      IF(SCFTYP.EQ.GVB) LGRDAF=402
      IF(SCFTYP.EQ.RMC) LGRDAF=403
      IF(CITYP.EQ.GUGA) LGRDAF=404
      CALL DAREAD(IDAF,IODA,E,L3ORB,LGRDAF,0)
C
C     ----- MAP MCSCF LAGRANGIAN INTO NORB BY NORB FORMAT -----
C
      IF(SCFTYP.EQ.GVB .OR. CITYP.EQ.GUGA) GO TO 280
C
      IF(NORB.GT.1) THEN
         IADD = 1
         IADDE = 1
         DO 260 I = 1,NORB
            IADD = IADD + NORB
            IADDE = IADDE + L1
            CALL DCOPY(NORB,E(IADDE),1,E(IADD),1)
  260    CONTINUE
      END IF
C
C     ----- CALCULATE C-DAGGER*E*C -----
C
  280 CONTINUE
      CALL VCLR(EPS,1,(L1*L1+L1)/2)
C
      IF (CITYP.NE.GUGA) THEN
C
C     ---CALCULATE THE HALF TRANSFORM FIRST -----
C     NOTE THAT E(KL) IS USED EXACTLY TWICE.  DIVIDE BY TWO TO GET
C     THE VALUES APPROPRIATE FOR THE GENERALIZED LAGRANGE MULTIPLIERS.
C
      DO 400 I = 1,L1
         KL = 0
         DO 330 L = 1,NORB
            WRK(L) = ZERO
            DO 320 K = 1,NORB
               KL = KL + 1
               WRK(L) = WRK(L) - V(I,K)*E(KL)
  320       CONTINUE
  330    CONTINUE
         DO 340 L = 1,NORB
            WRK(L) = WRK(L)*HALF
  340    CONTINUE
C
         DO 380 J = 1,L1
            IJ = IA(I) + J
            IF(J .GT. I) IJ = IA(J) + I
            DO 360 L = 1,NORB
               EPS(IJ) = EPS(IJ) + WRK(L)*V(J,L)
  360       CONTINUE
  380    CONTINUE
  400 CONTINUE
C
      ELSE
C
C     --- FOR CI THE LAGRANGIAN -ECI- IS NOT SYMMETRIC ---
C
      DO 500 I = 1,L1
         DO 430 L = 1,NORB
            WRK(L) = ZERO
            DO 420 K = 1,NORB
               IF(K.LE.L) THEN
                  WRK(L) = WRK(L) - V(I,K)*ECI(K,L)
               ELSE
                  WRK(L) = WRK(L) - V(I,K)*ECI(L,K)
               END IF
  420       CONTINUE
  430    CONTINUE
         DO 440 L = 1,NORB
            WRK(L) = WRK(L)*HALF
  440    CONTINUE
C
         DO 480 J = 1,L1
            IJ = IA(I) + J
            IF(J .GT. I) IJ = IA(J) + I
            DO 460 L = 1,NORB
               EPS(IJ) = EPS(IJ) + WRK(L)*V(J,L)
  460       CONTINUE
  480    CONTINUE
  500 CONTINUE
C
      END IF
C
C     CORRECT MCSCF AND GVB LAGRANGIANS. MP2 AND CI LARGANGIANS ARE
C     LEFT AS IS, SINCE THE METHODS ARE NOT FULLY VARIATIONAL.
C        V IS USED AS TEMPORARY STORAGE
C
      IF(NFG.NE.0.AND.(SCFTYP.EQ.GVB.OR.SCFTYP.EQ.RMC)) THEN
         L2 = (L1*L1+L1)/2
         IF(IAND(MODGRD,1).EQ.0) CALL FMOLAG(EPS,V,WRK,L1,L2)
      ENDIF
C
      RETURN
      END
C*MODULE GRD1    *DECK HELFEY
      SUBROUTINE HELFEY(DENAB,HFINT,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL IANDJ,OUT,DBG,NORM,DOUBLE,GOPARR,DSKWRK,MASWRK,NXT,SVDSKW
C
C  MAGIC NUMBERS: S(1)+P(3)+D(6)+F(10)+G(15)=35
C                 BASIS FUNCTION TYPES=5
C                 MAXIMUM BASIS TYPE SQUARED = (15)**2 = 225
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MXAO=8192)
C
      DIMENSION DENAB(L2),HFINT(L2,3,NAT)
      DIMENSION IJX(35),IJY(35),IJZ(35),DIJ(225),
     *          XIN(5,5,5,2),YIN(5,5,5,2),ZIN(5,5,5,2)
C
      COMMON /DSTV  / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,
     *                XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,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 /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /FMCOM / X(1)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      PARAMETER (RLN10=2.30258D+00)
      PARAMETER (PI212=1.1283791670955D+00)
      PARAMETER (SQRT3=1.73205080756888D+00)
      PARAMETER (SQRT5=2.23606797749979D+00)
      PARAMETER (SQRT7=2.64575131106459D+00)
C
      DATA DEBUG/8HDEBUG   /, DBUGME/8HHELFEY  /, GRD1/8HGRD1    /,
     *     RESC/8HRESC    /
C
C  THE IJX, IJY, AND IJZ ARRAYS CONTAIN THE POWERS OF THE CARTESIAN
C  GAUSSIANS PLUS 1 IN EVERY PLACE.  THE REASON FOR THE PLUS ONE IS
C  THAT THEY ARE USED TO INDEX AN ARRAY RUNNING FROM 1-5 CORRESPONDING
C  TO 0-4 IN THE POWER.
C
      DATA IJX/ 1, 2, 1, 1, 3, 1, 1, 2, 2, 1,
     *          4, 1, 1, 3, 3, 2, 1, 2, 1, 2,
     *          5, 1, 1, 4, 4, 2, 1, 2, 1, 3,
     *          3, 1, 3, 2, 2/
      DATA IJY/ 1, 1, 2, 1, 1, 3, 1, 2, 1, 2,
     *          1, 4, 1, 2, 1, 3, 3, 1, 2, 2,
     *          1, 5, 1, 2, 1, 4, 4, 1, 2, 3,
     *          1, 3, 2, 3, 2/
      DATA IJZ/ 1, 1, 1, 2, 1, 1, 3, 1, 2, 2,
     *          1, 1, 4, 1, 2, 1, 2, 3, 3, 2,
     *          1, 1, 5, 1, 2, 1, 2, 4, 4, 1,
     *          3, 3, 2, 2, 3/
C
C     ----- HELMANN-FEYNMAN GRADIENT TERM -----
C     INTEGRAL TYPE IS <II/H'/JJ> = <II/V'/JJ>
C     RESC RUNS DO NOT CONTRACT THE INTEGRALS WITH DENSITY
C
      DBG = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME)  .AND.  MASWRK
      OUT = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME  .OR.
     *       EXETYP.EQ.GRD1   .OR.  NPRINT.EQ.-3)
      IF(DBG) WRITE(IW,9000)
C
      TOL = RLN10*ITOL
      NORM = NORMF.NE.1 .OR. NORMP.NE.1
      IF(RMETHOD.EQ.RESC) CALL VCLR(HFINT,1,L2*3*NAT)
C
C     INITIALIZATION FOR PARALLEL
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT = -1
      MINE = -1
C
C     ----- I SHELL
C
      DO 600 II = 1,NSHELL
C
C           GO PARALLEL!
C
      IF(NXT .AND. GOPARR) THEN
         MINE = MINE + 1
         IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
         IF(NEXT.NE.MINE) GO TO 600
      END IF
C
      I = KATOM(II)
      XI = C(1,I)
      YI = C(2,I)
      ZI = C(3,I)
      I1 = KSTART(II)
      I2 = I1+KNG(II)-1
      LIT = KTYPE(II)
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II)-MINI
C
C     ----- J SHELL
C
      DO 580 JJ = 1,II
C
C           GO PARALLEL!
C
        IF((.NOT.NXT) .AND. GOPARR) THEN
           IPCOUNT = IPCOUNT + 1
           IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 580
        END IF
C
        J = KATOM(JJ)
        XJ = C(1,J)
        YJ = C(2,J)
        ZJ = C(3,J)
        J1 = KSTART(JJ)
        J2 = J1+KNG(JJ)-1
        LJT = KTYPE(JJ)
        MINJ = KMIN(JJ)
        MAXJ = KMAX(JJ)
        LOCJ = KLOC(JJ)-MINJ
        NROOTS = (LIT+LJT+1-2)/2 + 1
        RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
        IANDJ = II .EQ. JJ
C
C     ----- I PRIMITIVE
C
        DO 520 IG = I1,I2
        AI = EX(IG)
        ARRI = AI*RR
        AXI = AI*XI
        AYI = AI*YI
        AZI = AI*ZI
        CSI = CS(IG)
        CPI = CP(IG)
        CDI = CD(IG)
        CFI = CF(IG)
        CGI = CG(IG)
C
C     ----- J PRIMITIVE
C
        JGMAX = J2
        IF(IANDJ) JGMAX = IG
        DO 500 JG = J1,JGMAX
          AJ = EX(JG)
          AA = AI+AJ
          AA1 = ONE/AA
          DUM = AJ*ARRI*AA1
          IF(DUM .GT. TOL) GO TO 500
          FAC = EXP(-DUM)
          CSJ = CS(JG)
          CPJ = CP(JG)
          CDJ = CD(JG)
          CFJ = CF(JG)
          CGJ = CG(JG)
          AX = (AXI+AJ*XJ)*AA1
          AY = (AYI+AJ*YJ)*AA1
          AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
          DOUBLE=IANDJ.AND.IG.NE.JG
          JMAX = MAXJ
          NN = 0
          DUM1 = ZERO
          DUM2 = DUM1
          DO 360 I = MINI,MAXI
            IF(I.EQ.1) DUM1=CSI*FAC
            IF(I.EQ.2) DUM1=CPI*FAC
            IF(I.EQ.5) DUM1=CDI*FAC
            IF(I.EQ.8.AND.NORM) DUM1=DUM1*SQRT3
            IF(I.EQ.11) DUM1=CFI*FAC
            IF(I.EQ.14.AND.NORM) DUM1=DUM1*SQRT5
            IF(I.EQ.20.AND.NORM) DUM1=DUM1*SQRT3
            IF(I.EQ.21) DUM1=CGI*FAC
            IF(I.EQ.24.AND.NORM) DUM1=DUM1*SQRT7
            IF(I.EQ.30.AND.NORM) DUM1=DUM1*SQRT5/SQRT3
            IF(I.EQ.33.AND.NORM) DUM1=DUM1*SQRT3
C
            IF(IANDJ) JMAX = I
            DO 360 J = MINJ,JMAX
              IF(J.EQ.1) THEN
                DUM2=DUM1*CSJ
                IF( .NOT. DOUBLE) GO TO 350
                IF(I .GT. 1) THEN
                 DUM2 = DUM2+CSI*CPJ*FAC
                ELSE
                 DUM2 = DUM2+DUM2
                END IF
              ELSE IF(J.EQ.2) THEN
                DUM2=DUM1*CPJ
                IF(DOUBLE) DUM2 = DUM2+DUM2
              ELSE IF(J.EQ.5) THEN
                DUM2=DUM1*CDJ
                IF(DOUBLE) DUM2 = DUM2+DUM2
              ELSE IF(J.EQ.8.AND.NORM) THEN
                DUM2 = DUM2*SQRT3
              ELSE IF(J.EQ.11) THEN
                DUM2=DUM1*CFJ
                IF(DOUBLE) DUM2=DUM2+DUM2
              ELSE IF(J.EQ.14.AND.NORM) THEN
                DUM2=DUM2*SQRT5
              ELSE IF(J.EQ.20.AND.NORM) THEN
                DUM2=DUM2*SQRT3
              ELSE IF(J.EQ.21) THEN
                DUM2=DUM1*CGJ
                IF(DOUBLE) DUM2=DUM2+DUM2
              ELSE IF(J.EQ.24.AND.NORM) THEN
                DUM2=DUM2*SQRT7
              ELSE IF(J.EQ.30.AND.NORM) THEN
                DUM2=DUM2*SQRT5/SQRT3
              ELSE IF(J.EQ.33.AND.NORM) THEN
                DUM2=DUM2*SQRT3
              END IF
C
  350         NN = NN+1
              IF(RMETHOD.EQ.RESC) THEN
                DIJ(NN)=DUM2*PI212*AA1
              ELSE
                NDUM = IA(LOCI+I)+(LOCJ+J)
                DEN = DENAB(NDUM)
                IF(.NOT.IANDJ.OR.I.NE.J) DEN=DEN+DEN
                DIJ(NN)=DUM2*DEN*PI212*AA1
              ENDIF
  360     CONTINUE
C
C     ..... HELLMANN-FEYNMAN TERM .....
C
          AAX = AA*AX
          AAY = AA*AY
          AAZ = AA*AZ
          DO 480 IC = 1,NAT
            ZNUC = -ZAN(IC)
            CX = C(1,IC)
            CY = C(2,IC)
            CZ = C(3,IC)
            XX = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
            IF(NROOTS.LE.3) CALL RT123
            IF(NROOTS.EQ.4) CALL ROOT4
            IF(NROOTS.EQ.5) CALL ROOT5
            DO 420 K = 1,NROOTS
              UU = AA*U(K)
              WW = W(K)*ZNUC
              WW=WW*(UU+UU)
              TT = ONE/(AA+UU)
              T = SQRT(TT)
              X0 = (AAX+UU*CX)*TT
              Y0 = (AAY+UU*CY)*TT
              Z0 = (AAZ+UU*CZ)*TT
              DO 400 J = 1,LJT
                NJ = J
                DO 400 I = 1,LIT
                  NI = I
                  CALL VINT
                  XIN(I,J,K,1) = XINT
                  YIN(I,J,K,1) = YINT
                  ZIN(I,J,K,1) = ZINT*WW
                  CALL DVINT
                  XIN(I,J,K,2) = XINT
                  YIN(I,J,K,2) = YINT
                  ZIN(I,J,K,2) = ZINT*WW
  400         CONTINUE
  420       CONTINUE
           IJ=0
           DO 460 I=MINI,MAXI
            IX=IJX(I)
            IY=IJY(I)
            IZ=IJZ(I)
            JMAX=MAXJ
            IF(IANDJ) JMAX=I
            DO 460 J=MINJ,JMAX
             JX=IJX(J)
             JY=IJY(J)
             JZ=IJZ(J)
             DUMX = ZERO
             DUMY = ZERO
             DUMZ = ZERO
             DO 440 K = 1,NROOTS
              DUMX = DUMX+XIN(IX,JX,K,2)*YIN(IY,JY,K,1)*ZIN(IZ,JZ,K,1)
              DUMY = DUMY+XIN(IX,JX,K,1)*YIN(IY,JY,K,2)*ZIN(IZ,JZ,K,1)
  440         DUMZ = DUMZ+XIN(IX,JX,K,1)*YIN(IY,JY,K,1)*ZIN(IZ,JZ,K,2)
             IJ=IJ+1
             DUM=DIJ(IJ)
             IF(RMETHOD.EQ.RESC) THEN
               NN=IA(LOCI+I)+(LOCJ+J)
               HFINT(NN,1,IC)=HFINT(NN,1,IC)+DUMX*DUM
               HFINT(NN,2,IC)=HFINT(NN,2,IC)+DUMY*DUM
               HFINT(NN,3,IC)=HFINT(NN,3,IC)+DUMZ*DUM
             ELSE
               DE(1,IC) = DE(1,IC)+DUM*DUMX
               DE(2,IC) = DE(2,IC)+DUM*DUMY
               DE(3,IC) = DE(3,IC)+DUM*DUMZ
             ENDIF
             IF(DBG) WRITE(IW,9100) IC,II,JJ,IJ,DIJ(IJ),DUMX,DUMY,
     +                         DUMZ,DE(1,IC),DE(2,IC),DE(3,IC)
  460       CONTINUE
  480     CONTINUE
  500 CONTINUE
  520 CONTINUE
C
C     ----- END OF *PRIMITIVE* LOOPS -----
C
      IF(DBG) THEN
         WRITE(IW,9010) II,JJ
         CALL EGOUT(DE,NAT)
      END IF
  580 CONTINUE
  600 CONTINUE
C
C     ----- END OF *SHELL* LOOPS -----
C
      IF(OUT) THEN
         IF(GOPARR) CALL DDI_GSUMF(1502,DE,3*NAT)
         IF(MASWRK) THEN
            WRITE(IW,9000)
            CALL EGOUT(DE,NAT)
         END IF
         IF(GOPARR) CALL DSCAL(3*NAT,ONE/NPROC,DE,1)
      END IF
      IF(GOPARR .AND. NXT) CALL DDI_DLBRESET
      IF(RMETHOD.EQ.RESC) THEN
        IF(GOPARR) CALL DDI_GSUMF(1502,HFINT,L2*3*NAT)
        SVDSKW=DSKWRK
        DSKWRK=.TRUE.
        LREC=3*3*NAT
C       STARTING RECORD FOR THE HELLMANN-FEYNMAN INTEGRALS
        DO I=1,NAT
          DO J=1,3
            LREC=LREC+1
            CALL RAWRIT(NQRDAF,X(MORDA),HFINT(1,J,I),L2,LREC,0)
          ENDDO
        ENDDO
        DSKWRK=SVDSKW
      ENDIF
C
      RETURN
C
 9000 FORMAT(/10X,22(1H-)/10X,'HELLMANN-FEYNMAN FORCE'/10X,22(1H-))
 9010 FORMAT(1X,'HELFEY: SHELLS II,JJ=',2I5)
 9100 FORMAT(1X,'IC',I3,' II, JJ, I',3I5,1P,4E15.7/
     *       25X,'DE(1,2,3 AND IC)',1P,3E15.7)
      END
C*MODULE GRD1    *DECK SDER
      SUBROUTINE SDER(EPS,DSINT,L1,L2,NOCON)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DBG,NORM,OUT
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT,SVDSKW,NOCON
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MXAO=8192)
C
      DIMENSION EPS(L2),DSINT(L2,3,NAT)
      DIMENSION DIJ(225), DIJINT(225), IJX(35), IJY(35),IJZ(35),
     *           XS(6,5), YS(6,5), ZS(6,5),
     *          DXS(5,5),DYS(5,5),DZS(5,5)
C
      COMMON /DSTV  / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,
     *                XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,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 /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /FMCOM / X(1)
C
      PARAMETER (SQRT3=1.73205080756888D+00)
      PARAMETER (SQRT5=2.23606797749979D+00, SQRT7=2.64575131106459D+00)
      PARAMETER (ONE=1.0D+00)
      PARAMETER (RLN10=2.30258D+00)
C
      DATA DEBUG/8HDEBUG   /, DBUGME/8HSDER    /, GRD1/8HGRD1    /,
     *     RESC/8HRESC    /
C
      DATA IJX / 1, 2, 1, 1, 3, 1, 1, 2, 2, 1,
     *           4, 1, 1, 3, 3, 2, 1, 2, 1, 2,
     *           5, 1, 1, 4, 4, 2, 1, 2, 1, 3,
     *           3, 1, 3, 2, 2/
      DATA IJY / 1, 1, 2, 1, 1, 3, 1, 2, 1, 2,
     *           1, 4, 1, 2, 1, 3, 3, 1, 2, 2,
     *           1, 5, 1, 2, 1, 4, 4, 1, 2, 3,
     *           1, 3, 2, 3, 2/
      DATA IJZ / 1, 1, 1, 2, 1, 1, 3, 1, 2, 2,
     *           1, 1, 4, 1, 2, 1, 2, 3, 3, 2,
     *           1, 1, 5, 1, 2, 1, 2, 4, 4, 1,
     *           3, 3, 2, 2, 3/
C
C     ----- CALCULATE DERIVATIVES OF THE OVERLAP MATRIX -----
C     THIS TERM RESULTS FROM DERIVATIVES OF THE DENSITY MATRIX,
C     WHICH CAN IN THE CASE OF A SYMMETRIC LAGRANGIAN MATRIX, CAN
C     BE CONVERTED INSTEAD INTO DERIVATIVES OF OVERLAP INTEGRALS.
C     RESC RUNS NEED BOTH CONTRACTION OF THE INTEGRALS WITH DENSITY
C     AND THE INTEGRALS THEMSELVES.
C     IN CASE OF INTERNALLY UNCONTRACTED RESC, ONE HAS TO TRANSFORM
C     THE OVERLAP DERIVATIVES TO CONTRACTED BASIS BEFORE
C     SUMMING (CONTRACTING) OVERLAP DERIVATIVES WITH THE LAGRANGIAN
C
C     NOCON IS TRUE IF LAGRANGIAN IS NOT TO BE SUMMED WITH.
C
      DBG = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME)  .AND.  MASWRK
      OUT = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME  .OR.
     *       EXETYP.EQ.GRD1   .OR.  NPRINT.EQ.-3)
      IF(DBG) WRITE(IW,9000)
C
      TOL = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
      IF(RMETHOD.EQ.RESC) CALL VCLR(DSINT,1,L2*3*NAT)
C
C     INITIALIZE PARALLEL
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT = -1
      MINE = -1
C
      IF(DBG) THEN
         WRITE(IW,9010)
         CALL PRTRI(EPS,L1)
      END IF
C
C     ----- I SHELL
C
      DO 780 II = 1,NSHELL
C
C           GO PARALLEL!
C
      IF(NXT .AND. GOPARR) THEN
         MINE = MINE + 1
         IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
         IF(NEXT.NE.MINE) GO TO 780
      END IF
C
      IAT = KATOM(II)
      XI = C(1,IAT)
      YI = C(2,IAT)
      ZI = C(3,IAT)
      I1 = KSTART(II)
      I2 = I1+KNG(II)-1
      LIT = KTYPE(II)
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II)-MINI
      LITDER = LIT+1
C
C     ----- J SHELL
C
      DO 760 JJ = 1,II
C
C           GO PARALLEL!
C
      IF((.NOT.NXT) .AND. GOPARR) THEN
         IPCOUNT = IPCOUNT + 1
         IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 760
      END IF
C
      JAT = KATOM(JJ)
      XJ = C(1,JAT)
      YJ = C(2,JAT)
      ZJ = C(3,JAT)
      J1 = KSTART(JJ)
      J2 = J1+KNG(JJ)-1
      LJT = KTYPE(JJ)
      MINJ = KMIN(JJ)
      MAXJ = KMAX(JJ)
      LOCJ = KLOC(JJ)-MINJ
      RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
      IF(II.EQ.JJ) GO TO 760
C
C     ----- I PRIMITIVE
C
      DO 640 IG = I1,I2
        AI = EX(IG)
        ARRI = AI*RR
        AXI = AI*XI
        AYI = AI*YI
        AZI = AI*ZI
        CSI=CS(IG)
        CPI=CP(IG)
        CDI=CD(IG)
        CFI=CF(IG)
        CGI=CG(IG)
C
C     ----- J PRIMITIVE
C
        DO 620 JG = J1,J2
          AJ = EX(JG)
          AA = AI+AJ
          AA1 = ONE/AA
          DUM = AJ*ARRI*AA1
          IF(DUM .GT. TOL) GO TO 620
          FAC = EXP(-DUM)
          CSJ = CS(JG)
          CPJ = CP(JG)
          CDJ = CD(JG)
          CFJ = CF(JG)
          CGJ = CG(JG)
          AX = (AXI+AJ*XJ)*AA1
          AY = (AYI+AJ*YJ)*AA1
          AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
          IJ = 0
          DO 580 I = MINI,MAXI
            IF(I.EQ.1) DUM1=CSI*FAC
            IF(I.EQ.2) DUM1=CPI*FAC
            IF(I.EQ.5) DUM1=CDI*FAC
            IF(I.EQ.8.AND.NORM) DUM1=DUM1*SQRT3
            IF(I.EQ.11) DUM1=CFI*FAC
            IF(I.EQ.14.AND.NORM) DUM1=DUM1*SQRT5
            IF(I.EQ.20.AND.NORM) DUM1=DUM1*SQRT3
            IF(I.EQ.21) DUM1=CGI*FAC
            IF(I.EQ.24.AND.NORM) DUM1=DUM1*SQRT7
            IF(I.EQ.30.AND.NORM) DUM1=DUM1*SQRT5/SQRT3
            IF(I.EQ.33.AND.NORM) DUM1=DUM1*SQRT3
            DO 570 J = MINJ,MAXJ
              IF(J.EQ.1) DUM2=DUM1*CSJ
              IF(J.EQ.2) DUM2=DUM1*CPJ
              IF(J.EQ.5) DUM2=DUM1*CDJ
              IF(J.EQ.8.AND.NORM) DUM2=DUM2*SQRT3
              IF(J.EQ.11) DUM2=DUM1*CFJ
              IF(J.EQ.14.AND.NORM) DUM2=DUM2*SQRT5
              IF(J.EQ.20.AND.NORM) DUM2=DUM2*SQRT3
              IF(J.EQ.21) DUM2=DUM1*CGJ
              IF(J.EQ.24.AND.NORM) DUM2=DUM2*SQRT7
              IF(J.EQ.30.AND.NORM) DUM2=DUM2*SQRT5/SQRT3
              IF(J.EQ.33.AND.NORM) DUM2=DUM2*SQRT3
              IJ=IJ+1
              NN=IA(LOCI+I)+(LOCJ+J)
              IF(RMETHOD.EQ.RESC) THEN
                DIJINT(IJ)=DUM2
              ENDIF
              DEN = EPS(NN)
              DEN=DEN+DEN
              DIJ(IJ) = DUM2*DEN
  570       CONTINUE
  580     CONTINUE
C
C     ----- OVERLAP
C
      T = SQRT(AA1)
      X0 = AX
      Y0 = AY
      Z0 = AZ
      DO 590 J = 1,LJT
        NJ = J
        DO 590 I = 1,LITDER
          NI = I
          CALL VINT
          XS(I,J)=XINT*T
          YS(I,J)=YINT*T
          ZS(I,J)=ZINT*T
  590 CONTINUE
C
      CALL DERI(DXS,DYS,DZS,XS,YS,ZS,LIT,LJT,AI)
C
      IJ=0
      DO 600 I=MINI,MAXI
        IX=IJX(I)
        IY=IJY(I)
        IZ=IJZ(I)
        DO 600 J=MINJ,MAXJ
          JX=IJX(J)
          JY=IJY(J)
          JZ=IJZ(J)
          DUMX=DXS(IX,JX)* YS(IY,JY)* ZS(IZ,JZ)
          DUMY= XS(IX,JX)*DYS(IY,JY)* ZS(IZ,JZ)
          DUMZ= XS(IX,JX)* YS(IY,JY)*DZS(IZ,JZ)
          IJ=IJ+1
          IF(RMETHOD.EQ.RESC) THEN
            NN=IA(LOCI+I)+(LOCJ+J)
            DUMINT=DIJINT(IJ)
            DSINT(NN,1,IAT)=DSINT(NN,1,IAT)+DUMX*DUMINT
            DSINT(NN,2,IAT)=DSINT(NN,2,IAT)+DUMY*DUMINT
            DSINT(NN,3,IAT)=DSINT(NN,3,IAT)+DUMZ*DUMINT
            DSINT(NN,1,JAT)=DSINT(NN,1,JAT)-DUMX*DUMINT
            DSINT(NN,2,JAT)=DSINT(NN,2,JAT)-DUMY*DUMINT
            DSINT(NN,3,JAT)=DSINT(NN,3,JAT)-DUMZ*DUMINT
          ENDIF
          IF(.NOT.NOCON) THEN
            DE(1,IAT)=DE(1,IAT)+(DUMX*DIJ(IJ))
            DE(2,IAT)=DE(2,IAT)+(DUMY*DIJ(IJ))
            DE(3,IAT)=DE(3,IAT)+(DUMZ*DIJ(IJ))
            DE(1,JAT)=DE(1,JAT)-(DUMX*DIJ(IJ))
            DE(2,JAT)=DE(2,JAT)-(DUMY*DIJ(IJ))
            DE(3,JAT)=DE(3,JAT)-(DUMZ*DIJ(IJ))
          ENDIF
  600 CONTINUE
C
  620 CONTINUE
  640 CONTINUE
C
C     ----- END OF PRIMITIVE LOOPS -----
C
      IF(DBG) THEN
         WRITE(IW,9100) II,JJ
         CALL EGOUT(DE,NAT)
      END IF
  760 CONTINUE
  780 CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
      IF(OUT) THEN
         IF(GOPARR) CALL DDI_GSUMF(1501,DE,3*NAT)
         IF(MASWRK) THEN
            WRITE(IW,9000)
            CALL EGOUT(DE,NAT)
         END IF
         IF(GOPARR) CALL DSCAL(3*NAT,ONE/NPROC,DE,1)
      END IF
      IF(GOPARR .AND. NXT) CALL DDI_DLBRESET
C
      IF(RMETHOD.EQ.RESC) THEN
        IF(GOPARR) CALL DDI_GSUMF(1501,DSINT,L2*3*NAT)
        SVDSKW=DSKWRK
        DSKWRK=.TRUE.
        LREC=0
C       STARTING RECORD FOR THE S DERIVATIVE INTEGRALS
        DO I=1,NAT
          DO J=1,3
            LREC=LREC+1
            CALL RAWRIT(NQRDAF,X(MORDA),DSINT(1,J,I),L2,LREC,0)
          ENDDO
        ENDDO
        DSKWRK=SVDSKW
      ENDIF
C
      RETURN
C
 9000 FORMAT(/10X,33(1H-)/10X,'GRADIENT INCLUDING DENSITY FORCES'/
     *        10X,33(1H-))
 9010 FORMAT(1X,'THE LAGRANGIAN IN THE AO BASIS IS')
 9100 FORMAT(1X,'SDER: SHELLS II,JJ=',2I5)
      END
C*MODULE GRD1    *DECK STVDER
      SUBROUTINE STVDER
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL SOME,OUT,GOPARR,DSKWRK,MASWRK,EFLDL,SVDSKW,ISEPS,USEPS,
     *        UNCON,NXT
C
      PARAMETER (MXATM=2000, MXAO=8192, MXSH=5000, MXGTOT=20000,
     *           MXPT=2000, MXFRG=1050, MXFGPT=12000, MXDFG=5,
     *           MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (MAXDEN=25*MXATM, LENABC=2000)
C
      COMMON /EFLDC / EVEC(3),EFLDL
      COMMON /EFPBAS/ EX(MXGTOT,MXDFG),CS(MXGTOT,MXDFG),
     *                CP(MXGTOT,MXDFG),CD(MXGTOT,MXDFG),
     *                CF(MXGTOT,MXDFG),CG(MXGTOT,MXDFG),
     *                PRNAME(MXFGPT),PRCORD(3,MXFGPT),EFZNUC(MXFGPT),
     *                KSTART(MXSH,MXDFG),KATOM(MXSH,MXDFG),
     *                KTYPE(MXSH,MXDFG),KNG(MXSH,MXDFG),
     *                KLOC(MXSH,MXDFG),KMIN(MXSH,MXDFG),
     *                KMAX(MXSH,MXDFG),NSHELL(MXDFG),NGAUSS(MXDFG),
     *                NATEF(MXDFG),NUMEF(MXDFG),NTPATM
      COMMON /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      COMMON /EFPOTD/ METHOF,ISCCHG,MOVE
      COMMON /ELGPMT/ NELONG,NATM,NASPIN,NCT,NBNDAB,NTMLB,IPRI,LDOS
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FMCOM / X(1)
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /ISEPS / ISEPS, USEPS
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB(MXFRG),
     *                NPBF(MXFRG),NTMO
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SOLV  / FEPSI,RDS,DISEX2,COSURF(3,LENABC),SRAD(MXATM),
     *                QDEN(MAXDEN),AR(LENABC),
     *                NSPA,NPS,NPS2,NDEN,NPSPHER
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
      DATA CHECK/8HCHECK   /, ZRF/8HZRF     /,GRD1/8HGRD1    /
      DATA DEBUG/8HDEBUG   /, ANONE/8HNONE    /
      DATA ANESC,RESC,DK/8HNESC    ,8HRESC    ,8HDK      /
C
C     ----- DRIVER FOR THE GRADIENT OF THE ONE ELECTRON TERMS -----
C
      IPOT=IEFC+IEFD+IEFQ+IEFO+IEFP+IEFDP+IREP+NTMO
C
      SOME = NPRINT.NE.-5 .AND. MASWRK
      OUT = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.GRD1   .OR.  NPRINT.EQ.-3)
     *         .AND.  MASWRK
C
C     THERE IS NO DK ANALYTIC GRADIENT,
C     SO GRADIENT DK CALCULATIONS DO NOT NEED UNCONTRACTION.
C
      UNCON=RMETHOD.EQ.RESC.AND.MOD(MODQR,2).EQ.1
C
      IF(SOME) THEN
         T0 = 0.0D+00
         IF(IECP.GT.0  .OR.  IZRF.GT.0  .OR.  IPOT.GT.0 .OR. EFLDL)
     *     CALL TSECND(T0)
           WRITE(IW,9000)
      END IF
C
C     INITIALIZE GRADIENT TO ZERO
C
      NFRPTS = NMTTPT + NPTTPT + NRTTPT + NTPATM + NTMO + NDTTPT
      CALL VCLR(DE ,1,3*NAT)
      CALL VCLR(DEF,1,3*MXFGPT)
      IF(NFRG.GT.0) CALL VCLR(ATORQ,1,3*NFRG)
      IF(NTMO.GT.0) CALL DPAULA
C
C     ----- GET WORKING STORAGE -----
C     WE REALLY WANT -LEIG- TO BE A SQUARE AS THAT IS WORKING STORE
C     FOR SOME OF THE LAGRANGIAN SETUP.
C
      L1 = NUM
      L2 = (NUM*NUM+NUM)/2
      L3 = NUM*NUM
C
      CALL VALFM(LOADFM)
      LEPS   = LOADFM + 1
      LVEC   = LEPS   + L2
      LEIG   = LVEC   + L3
      LIA    = LEIG   + L3
      LWRK   = LIA    + L1
      LDA    = LWRK   + L1
      LDB    = LDA    + L2
      LDRG   = LDB    + L2
      LAST   = LDRG   + NAT*NAT
C
      IF(NELONG.GT.2) THEN
        LFA  = LAST
        LFX  = LFA    + L2
        LAST = LFX    + L3
      ELSE
        LFA  = LAST
        LFX  = LAST
      ENDIF
C
C          EXTRA STORAGE FOR SOME SOLVATION MODELS
                NDIMCD=1
      IF(ISEPS) NDIMCD=NAT+NPS
      LEF3   = LAST
      LCOSDG = LEF3   + MAX(1,3*NFRPTS)
      LWORK  = LCOSDG + NDIMCD*NDIMCD
      LAST   = LWORK  + MAX(1,10*NMTTPT)
C
                    MXTMP = 1
      IF(NFRG.GT.0) MXTMP = MXFGPT
      LEFATMP = LAST    + 3*MXTMP
      LEFBTMP = LEFATMP + 3*MXTMP
      LCHGTMP = LEFBTMP + 3*MXTMP
      LDIPTMP = LCHGTMP + 3*MXTMP
      LQUATMP = LDIPTMP + 3*MXTMP
      LINDTMP = LQUATMP + 3*MXTMP
      LIDDTMP = LINDTMP + 3*MXTMP
      LEF3FRG = LIDDTMP + 3*MXTMP
      LATORQN = LEF3FRG + 3*MXTMP
      LAST    = LATORQN + MAX(1,3*NFRG)
C
C          EXTRA STORAGE FOR RELATIVITY
      IF(RMETHOD.EQ.RESC  .OR.  RMETHOD.EQ.DK) THEN
         L2U=L2
         IF(UNCON) L2U=(NUMU*NUMU+NUMU)/2
         NQRDAF=30
         NTOTREC=4*3*NAT
         LDERIV = LAST
         LDERIV1= LDERIV  + L2U*3*NAT
         MORDA  = LDERIV1 + L2U*3*NAT
         LWRKL2 = MORDA   + (NTOTREC-1)/NWDVAR+1
         LAST   = LWRKL2  + L2
         IF(UNCON) THEN
            LUU  = LAST
            LWRKU= LUU   + NUMU*L1
            LAST = LWRKU + NUMU
         ELSE
            LUU  = LAST
            LWRKU= LAST
         ENDIF
      ELSE
         LDERIV =LAST
         LDERIV1=LAST
         MORDA  =LAST
         LWRKL2 =LAST
         LUU    =LAST
         LWRKU  =LAST
      ENDIF
C
C          EXTRA STORAGE FOR MCP DERIVATIVES
C              THE 900'S ARE FIRMLY EMBEDDED IN THE MCP CODE, SO YOU'D
C              NEED TO DO MORE THAN JUST CHANGE THE TWO SIZES HERE.
      MXVO  =900
      MXCO  =900
      IF(IECP.EQ.5) THEN
         LDD    = LAST
         LOV    = LDD  + 3*NAT
         LDV    = LOV  + MXCO*MXVO
         LAST   = LDV  + MXCO*MXVO*3
      ELSE
         LDD = LAST
         LOV = LAST
         LDV = LAST
      END IF
C
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      IF(EXETYP.EQ.CHECK) GO TO 100
C
C     ----- GET LAGRANGIAN -EPS- MATRIX -----
C
      IF(NELONG.LE.2) THEN
         CALL EIJDEN(X(LEPS),X(LVEC),X(LEIG),
     *               X(LIA),X(LWRK),L1,L2,L3,NQMT,X(LEIG))
      ELSE
         CALL ELGLAG(X(LEPS),X(LVEC),X(LEIG),X(LFA),X(LFX),L1,L2,NE)
      END IF
      IF(TDDFTYP.NE.ANONE) CALL TDDENS(X(LEPS),X(LVEC),L1)
      IF(NTMO.GT.0) CALL QMADD(X(LEPS),X(LDA),L2,78)
C
      IF(OUT) THEN
         WRITE(IW,*) 'THE LAGRANGIAN MATRIX IS'
         CALL PRTRI(X(LEPS),L1)
      END IF
C
C     ----- GET TOTAL DENSITY -DA- MATRIX -----
C
      CALL DENDD1(X(LDA),X(LDB),L2)
      IF(OUT) THEN
         WRITE(IW,*) 'THE DENSITY MATRIX IS'
         CALL PRTRI(X(LDA),L1)
      END IF
C
      IF(RMETHOD.EQ.RESC) THEN
         SVDSKW=DSKWRK
         DSKWRK=.TRUE.
         CALL RAOPEN(NQRDAF,X(MORDA),0,NTOTREC,L3,NPRINT)
         DSKWRK=SVDSKW
      ENDIF
C
      IF(RMETHOD.EQ.ANESC) CALL FLIPBASIS(0)
      IF(UNCON) THEN
         CALL FLIPBASIS(17)
         L1=NUMU
         L2=(L1*L1+L1)/2
         L3=L1*L1
      ENDIF
C
C     ----- NUCLEAR REPULSION FORCE -----
C     ----- DENSITY FORCE -----
C
      CALL VNNDER(X(LDRG),NAT,X(LEF3),NFRPTS,X(LCOSDG),NDIMCD,X(LWORK),
     *            X(LEFATMP),X(LEFBTMP),X(LCHGTMP),X(LDIPTMP),
     *            X(LQUATMP),X(LINDTMP),X(LIDDTMP),X(LEF3FRG),
     *            X(LATORQN),MXTMP)
      CALL SDER(X(LEPS),X(LDERIV),L1,L2,UNCON)
      IF (NTMO.GT.0) CALL QMADD(X(LDA),X(LEPS),L2,79)
C
C       QMADD CAN FAIL TO WORK WITH UNCONTRACTED BASIS (?)
C       IN THIS CASE SDER COULD NOT CONTRACT LAGRANGIAN SINCE
C       IT IS COMPUTED IN THE UNCONTRACTED SPACE AND THE DENSITY
C       IN CONTRACTED.
C       READ IN THE TRANSFORMATION MATRIX (CONTRACTION COEFFICIENTS)
C
      IF(UNCON) THEN
        CALL DAREAD(IDAF,IODA,X(LUU),L1*NUM,NDARELB+19,0)
        LI=LDERIV
        NXT = IBTYP.EQ.1
        NEXT = -1
        MINE = -1
        DO IATM=1,NAT
          DO IXYZ=1,3
            MINE = MINE + 1
            IF(GOPARR.AND.NXT) THEN
              IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
              IF(NEXT.NE.MINE) GO TO 50
            END IF
            IF(GOPARR.AND..NOT.NXT.AND.MOD(MINE,NPROC).NE.ME) GOTO 50
            CALL GRADEN(DE(IXYZ,IATM),X(LI),X(LWRKL2),X(LUU),X(LEPS),
     *                  X(LWRKU),UNCON,0)
   50       CONTINUE
            LI=LI+L2
          ENDDO
        ENDDO
        IF(GOPARR.AND.NXT) CALL DDI_DLBRESET
      ENDIF
C
C     ----- HELLMANN-FEYNMAN FORCE -----
C     ----- INTEGRAL FORCE (AO DERIVATIVE CONTRIBUTION) -----
C
      CALL HELFEY(X(LDA),X(LDERIV),L2)
      CALL TVDER(X(LDA),X(LDERIV),X(LDERIV1),L2,0)
C
C     DERIVATIVES OF PVP INTEGRALS (HELLMANN-FEYNMAN AND AO DERIVATIVE)
C     CALLS MUST BE IN THIS ORDER SINCE THE ARRAY IS ACCUMULATED
C
      IF(RMETHOD.NE.ANONE) THEN
         IF(RMETHOD.EQ.ANESC) CALL FLIPBASIS(1)
         IF(RMETHOD.EQ.ANESC  .OR.  RMETHOD.EQ.RESC) THEN
            CALL HFPVP(X(LDA),X(LDERIV),L2)
            CALL DPVPINT(X(LDA),X(LDERIV),L2)
         END IF
C
C        NESC ADDS CORRECTIONS TO S AND T
C
         IF(RMETHOD.EQ.ANESC) THEN
C                GET CORRECTION TO S (USE LAGRANGIAN!)
            CALL TVDER(X(LEPS),X(LDERIV),X(LDERIV1),L2,1)
            CALL FLIPBASIS(2)
C                GET CORRECTION TO T (USE DENSITY!)
            CALL TVDER(X(LDA),X(LDERIV),X(LDERIV1),L2,2)
            CALL FLIPBASIS(3)
         ENDIF
C
         IF(RMETHOD.EQ.RESC) CALL RESCGX(X(LDA),L1,L2,L3)
         IF(RMETHOD.EQ.DK) CALL QRNUMDER(X(LDA),L2)
C
C        THERE IS NO GRADIENT DRIVER FOR NESC - IT IS NOT NEEDED!
C
         IF(UNCON) THEN
            CALL FLIPBASIS(15)
            L1=NUM
            L2 = (L1*L1+L1)/2
            L3=L1*L1
         ENDIF
         IF(SOME.AND.(IECP.GT.0.OR.IZRF.GT.0.OR.IPOT.GT.0.OR.EFLDL))THEN
           CALL TSECND(T1)
           TR = T1-T0
           WRITE(IW,9015) TR
           T0 = T1
         END IF
      ENDIF
      CALL DENDD1(X(LDA),X(LDB),L2)
C
      IF(SOME  .AND.
     *  (IECP.GT.0 .OR. IZRF.GT.0 .OR. IPOT.GT.0 .OR. EFLDL)) THEN
         CALL TSECND(T1)
         TG = T1-T0
         WRITE(IW,9020) TG
         T0 = T1
      END IF
C
C     ----- EFFECTIVE CORE POTENTIAL MODIFICATIONS -----
C
C     *** MODEL CORE POTENTIAL POTENTIAL MODIFICATIONS WERE
C     *** ADDED BY HIROTOSHI MORI & EISAKU MIYOSHI
C
  100 CONTINUE
      IF(IECP.GT.0) THEN
         IF (IECP.EQ.5) THEN
            CALL MCPGRD(DE, X(LDA), X(LDD), X(LOV), X(LDV), MXCO, MXVO)
         ELSE
            CALL ECPGRD(DE, X(LDA), EXETYP, SOME)
         ENDIF
         CALL TSECND(T0)
      END IF
C
C     ----- GET SELF-CONSISTENT REACTION FIELD GRADIENT -----
C
      IF(IZRF.EQ.1) THEN
         DUMMY = 0.0D+00
         CALL DIPDER(X(LDA),DUMMY,DUMMY,DUMMY,L2,DUMMY,0,0,0,ZRF)
         IF(SOME) THEN
            CALL TSECND(T1)
            TG = T1-T0
            WRITE(IW,9040) TG
            T0 = T1
         END IF
      END IF
C
C     ----- EFFECTIVE FRAGMENT GRADIENT -----
C     EFFT CALCULATES TOTAL FORCE AND TORQUE ON EACH FRAGMENT
C
      IF(IPOT.GT.0) THEN
C
C   IF NBUFMO.GT.0 THEN WE HAVE TO USE THE BUFFERLESS DENSITY,
C   WHICH REQUIRES MORE MEMORY
C
         IF(NBUFMO.GT.0)THEN
            LVF   = LAST
            LDBF  = LVF   + L3
            LOCCF = LDBF  +  L2
            LWRKF = LOCCF +  L1
            LASTF = LWRKF +  L2
            NEEDF = LASTF - LAST
            CALL GETFM(NEEDF)
            CALL DCOPY(L2,X(LDA),1,X(LWRKF),1)
            CALL DENDBF(X(LDA),X(LDB),X(LDBF),X(LVF),X(LOCCF),NBUFMO,
     *                L1,L2,L3)
         END IF
         CALL FRGDER(X(LDA),L2)
         IF(NBUFMO.GT.0) THEN
             CALL DCOPY(L2,X(LWRKF),1,X(LDA),1)
             CALL RETFM(NEEDF)
         END IF
         IF(MOVE.NE.0) CALL EFFT
         IF(SOME) THEN
            CALL TSECND(T1)
            TG = T1-T0
            WRITE(IW,9050) TG
            T0 = T1
         END IF
      END IF
C
C     ----- GET ELECTRIC FIELD CONTRIBUTION -----
C
      IF (EFLDL) THEN
         CALL EFLD2(X(LDA),L2)
         IF(SOME) THEN
            CALL TSECND(T1)
            TG = T1-T0
            WRITE(IW,9060) TG
            T0 = T1
         END IF
      END IF
C
      IF(RMETHOD.EQ.RESC) CALL RACLOS(NQRDAF,'KEEP')
      CALL RETFM(NEED)
C
C     ----- GATHER ALL GRADIENT CONTRIBUTIONS -----
C     ----- SAVE COMPLETED 1E-GRADIENT -----
C
      IF(GOPARR) CALL DDI_GSUMF(1500,DE,3*NAT)
      CALL DAWRIT(IDAF,IODA,DE,3*NAT,3,0)
C
C     ----- CHECK TIME -----
C
      IREST=4
      IST=1
      JST=1
      KST=1
      LST=1
      IF(NFG.EQ.0.OR.NPRINT.NE.-5) THEN
        IF(MASWRK) WRITE(IW,9010)
        CALL TEXIT(1,IREST)
      ENDIF
      RETURN
C
 9000 FORMAT(/1X,'BEGINNING ONE ELECTRON GRADIENT...')
 9010 FORMAT(1X,'..... END OF 1-ELECTRON GRADIENT ......')
 9015 FORMAT(1X,'TIME TO DO RELATIVISTIC GRADIENT INTEGRALS=',F10.2)
 9020 FORMAT(1X,'TIME TO DO ORDINARY GRADIENT INTEGRALS=',F10.2)
 9040 FORMAT(1X,'TIME TO DO     SCRF GRADIENT INTEGRALS=',F10.2)
 9050 FORMAT(1X,'TIME TO DO FRAGMENT GRADIENT INTEGRALS=',F10.2)
 9060 FORMAT(1X,'TIME TO DO ELECTRIC FIELD CONTRIBUTION=',F10.2)
      END
C*MODULE GRD1    *DECK TVDER
      SUBROUTINE TVDER(DENAB,DTINT,DVINT,L2,MODUS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL OUT,DBG,NORM
      LOGICAL GOPARR,DSKWRK,MASWRK,NXT,SVDSKW,ESD
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MXAO=8192)
      PARAMETER (MXCHRM=1)
C
      DIMENSION DENAB(L2),DTINT(L2,3,NAT),DVINT(L2,3,NAT)
      DIMENSION DIJ(225), IJX(35), IJY(35), IJZ(35),
     *          XS(6,7), YS(6,7), ZS(6,7), XT(6,5), YT(6,5), ZT(6,5),
     *          DXS(5,5), DYS(5,5), DZS(5,5),
     *          DXT(5,5), DYT(5,5), DZT(5,5),
     *          XV(6,5,5), YV(6,5,5), ZV(6,5,5),
     *          DXV(5,5,5), DYV(5,5,5), DZV(5,5,5)
C
      COMMON /CHMGMS/ XCHM(MXCHRM),YCHM(MXCHRM),ZCHM(MXCHRM),
     *                DXELMM(MXCHRM),DYELMM(MXCHRM),DZELMM(MXCHRM),
     *                QCHM(MXCHRM),NCHMAT,KCHRMM
      COMMON /DSTV  / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,
     *                XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      COMMON /FMCOM / X(1)
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /FMOOPT/ ESPSCA(9),RESPAP(2),RESPPC(2),RESDIM,RESTRI(4),
     *                RCORSD,RESPCT,CONVFG,CNVDMP,COROFF,IXESP,MXITFG,
     *                NGUESS,NBSSE,MODORB,MODPAR,IRSTSTP,IRSTLAY,NPRFMO,
     *                NFMOPAL,MODPRP,MAXL1C,IPIEDA,MODGRD,MODESP,IVMUL
      COMMON /FMOPNT/ LICHFG,LMULFG,LIDMREC,LFRGNAM,LLAYFRG,LINDAT,
     *                LNCBS,LFMOZAN,LFMOC,LFMOMAS,LIZBAS,LIAGLOB,LIBDGH,
     *                LIABDFG,LJABDFG,LNCAO,LIDXCAO,LIAPRJO,LJAPRJO,
     *                LCOREAO,LOCCCOR,LSHIFTB,LIODFMO,LFMODA,LFMODB,
     *                LFMOESPA,LFMOESPB,LLOCFMO,LSCFFRG,LRIJ,LPOPMUL,
     *                LPOPMAT,LIALOC,LINDBD,LIATFRG,LINDFRG,LINDGFRG,
     *                LNATFRG,LNAT0FRG,LIANFRG,LZANFRG,LCFRG,LLIBISH,
     *                LLIBNSH,LLIBNG,LINDATG,LFMOBUF(3),LFMODE,LNUMFRG,
     *                LLOCTAT,LIAOGLOB,LLOADM,LFMOGE,LDGRID,LIODCFMO,
     *                LJOB2GRP,LFMOPG,LEMOCDR,LUNTXYZ,LUNTROT,LSTONEP,
     *                LMAPSU,LFRGMUL
      COMMON /FMORUN/ ESPSCF,E0SCF(2),IDAFMO,ICURLAY,ICURFG,JCURFG,
     *                KCURFG,NCURSH,NGAU,ICURPOP,IFMOSTP,MONCOR,NAT1E,
     *                NEEDR,MODRST,NORBPROJ,NUNESP,ISKIPESP,ICURUNT
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,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 /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (MAXDEN=25*MXATM, LENABC=2000, NPPA=1082)
      LOGICAL ISEPS, USEPS
      COMMON /ISEPS / ISEPS, USEPS
      COMMON /SOLV  / FEPSI,RDS,DISEX2,COSURF(3,LENABC),SRAD(MXATM),
     *                QDEN(MAXDEN),AR(LENABC),
     *                NSPA,NPS,NPS2,NDEN,NPSPHER
      COMMON /SOLVI / IATSP(LENABC+1),N0(2),NP1,NP2,ISKIP
      COMMON /COSMO3/ COSZAN(NPPA),CORZAN(3,NPPA)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00,TWO=2.0D+00)
      PARAMETER (RLN10=2.30258D+00)
      PARAMETER (PI212=1.1283791670955D+00)
      PARAMETER (SQRT3=1.73205080756888D+00)
      PARAMETER (SQRT5=2.23606797749979D+00,SQRT7=2.64575131106459D+00)
C
      DATA DEBUG/8HDEBUG   /, DBUGME/8HTVDER   /, GRD1/8HGRD1    /,
     *     RESC/8HRESC    /,ANESC/8HNESC    /
C
      DATA IJX / 1, 2, 1, 1, 3, 1, 1, 2, 2, 1,
     *           4, 1, 1, 3, 3, 2, 1, 2, 1, 2,
     *           5, 1, 1, 4, 4, 2, 1, 2, 1, 3,
     *           3, 1, 3, 2, 2/
      DATA IJY / 1, 1, 2, 1, 1, 3, 1, 2, 1, 2,
     *           1, 4, 1, 2, 1, 3, 3, 1, 2, 2,
     *           1, 5, 1, 2, 1, 4, 4, 1, 2, 3,
     *           1, 3, 2, 3, 2/
      DATA IJZ / 1, 1, 1, 2, 1, 1, 3, 1, 2, 2,
     *           1, 1, 4, 1, 2, 1, 2, 3, 3, 2,
     *           1, 1, 5, 1, 2, 1, 2, 4, 4, 1,
     *           3, 3, 2, 2, 3/
C
C     ----- BASIS FUNCTION DERIVATIVE CONTRIBUTIONS TO GRADIENT -----
C     INTEGRALS ARE OF TYPE <II'/H/JJ> = <II'/T+V/JJ>
C     RESC RUNS STORE INTEGRALS WITHOUT CONTRACTING THEM WITH DENSITY
C
C     MODUS - USED ONLY WITH NESC METHOD AND ONLY FOR KIN. ENERGY T
C           = 0 BULK T GRADIENT
C           = 1 CORRECTION TO S DERIVATIVE (PROPORTIONAL TO T)
C           = 2 CORRECTION TO T DERIVATIVE
C
      DBG = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME)  .AND.  MASWRK
      OUT = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME  .OR.
     *       EXETYP.EQ.GRD1   .OR.  NPRINT.EQ.-3)
      IF(DBG) WRITE(IW,9000)
      IAZ=0
C
      TOL = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
      IF(RMETHOD.EQ.ANESC) THEN
         FSC=ONE/CLIG
         FSC2=FSC*FSC/TWO
      ENDIF
      IF(RMETHOD.EQ.RESC) THEN
         CALL VCLR(DTINT,1,L2*3*NAT)
         CALL VCLR(DVINT,1,L2*3*NAT)
      ENDIF
C
C INITIALIZE PARALLEL
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT = -1
      MINE = -1
C
C     -D*ELMM- CONTAINS CONTRIBUTION TO THE FORCES ON THE
C     EXTERNAL CHARGES WHICH ARE DUE TO QM MOLECULE.
C
      IF(NCHMAT.NE.0) THEN
         DO 100 ICHARM=1,NCHMAT
            DXELMM(ICHARM)=ZERO
            DYELMM(ICHARM)=ZERO
            DZELMM(ICHARM)=ZERO
  100    CONTINUE
      END IF
      ESD=NFG.NE.0.AND. IFMOSTP.EQ.2.AND.MPLEVL.EQ.0.AND.
     *                  RESDIM.NE.0.AND.IAND(IXESP,32).EQ.0
      NEXTAT=NCHMAT
      IF(ESD) NEXTAT=NATFMO
C
C     FIRST, LOOP OVER THE GENUINE QM MOLECULE (ICHARM=0 PASS),
C     AND THEN ANY EXTERNAL CHARGES WHICH CHARMM MAY BE USING.
C     THE LATTER MAKE A CONTRIBUTION TO THE GRADIENT OF THE ATOMS
C     IN THE QM MOLECULE, AS WELL AS TO THEMSELVES.
C     FMO EXTERNAL ATOMS FROM SEPARATED DIMERS ARE HANDLED SIMILARLY TO
C     CHARMM.
C
      DO 2000 ICHARM=0,NEXTAT
        IF(ESD.AND.ICHARM.GT.0) THEN
          CALL ESDDERZ(ICHARM,X(LLAYFRG),X(LINDAT),X(LIABDFG),X(LJABDFG)
     *                ,X(LINDBD),X(LFMOZAN),X(LFMOC),ZNUC,CX,CY,CZ)
C         WRITE(6,*) 'WWWEXTZ',ICHARM,-ZNUC,CX,CY,CZ
          IF(ABS(ZNUC).LT.1.0D-08) GOTO 2000
C       THE CHARGE IS OF COURSE EITHER 0 OR NOT, NOT MERELY CLOSE TO 0.
      ENDIF
C
C     ----- I SHELL
C
      DO 1400 II = 1,NSHELL
C
C           GO PARALLEL!
C
        IF(NXT .AND. GOPARR) THEN
           MINE = MINE + 1
           IF(MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
           IF(NEXT.NE.MINE) GO TO 1400
        END IF
C
        IAT = KATOM(II)
        XI = C(1,IAT)
        YI = C(2,IAT)
        ZI = C(3,IAT)
        I1 = KSTART(II)
        I2 = I1+KNG(II)-1
        LIT = KTYPE(II)
        MINI = KMIN(II)
        MAXI = KMAX(II)
        LOCI = KLOC(II)-MINI
        LITDER = LIT + 1
C
C     ----- J SHELL
C
        DO 1300 JJ = 1,NSHELL
C
C           GO PARALLEL!
C
          IF((.NOT.NXT) .AND. GOPARR) THEN
             IPCOUNT = IPCOUNT + 1
             IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 1300
          END IF
C
          JAT = KATOM(JJ)
          XJ = C(1,JAT)
          YJ = C(2,JAT)
          ZJ = C(3,JAT)
          J1 = KSTART(JJ)
          J2 = J1+KNG(JJ)-1
          LJT = KTYPE(JJ)
          MINJ = KMIN(JJ)
          MAXJ = KMAX(JJ)
          LOCJ = KLOC(JJ)-MINJ
          LJTMOD = LJT+2
          NROOTS = (LIT+LJT-1)/2 + 1
          RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
C
C     ----- I PRIMITIVE
C
        DO 1100 IG = I1,I2
          AI = EX(IG)
          ARRI = AI*RR
          AXI = AI*XI
          AYI = AI*YI
          AZI = AI*ZI
          CSI=CS(IG)
          CPI=CP(IG)
          CDI=CD(IG)
          CFI=CF(IG)
          CGI=CG(IG)
C
C     ----- J PRIMITIVE
C
          DO 1000 JG = J1,J2
            AJ = EX(JG)
            AA = AI+AJ
            AA1 = ONE/AA
            DUM = AJ*ARRI*AA1
            IF(DUM .GT. TOL) GO TO 1000
            FAC = EXP(-DUM)
            CSJ = CS(JG)
            CPJ = CP(JG)
            CDJ = CD(JG)
            CFJ = CF(JG)
            CGJ = CG(JG)
            AX = (AXI+AJ*XJ)*AA1
            AY = (AYI+AJ*YJ)*AA1
            AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
            IJ = 0
            DUM1 = ZERO
            DUM2 = DUM1
            DO 600 I=MINI,MAXI
              IF(I.EQ.1) DUM1=CSI*FAC
              IF(I.EQ.2) DUM1=CPI*FAC
              IF(I.EQ.5) DUM1=CDI*FAC
              IF(I.EQ.8.AND.NORM) DUM1=DUM1*SQRT3
              IF(I.EQ.11) DUM1=CFI*FAC
              IF(I.EQ.14.AND.NORM) DUM1=DUM1*SQRT5
              IF(I.EQ.20.AND.NORM) DUM1=DUM1*SQRT3
              IF(I.EQ.21) DUM1=CGI*FAC
              IF(I.EQ.24.AND.NORM) DUM1=DUM1*SQRT7
              IF(I.EQ.30.AND.NORM) DUM1=DUM1*SQRT5/SQRT3
              IF(I.EQ.33.AND.NORM) DUM1=DUM1*SQRT3
C
              DO 600 J = MINJ,MAXJ
                IF(J.EQ.1) DUM2=DUM1*CSJ
                IF(J.EQ.2) DUM2=DUM1*CPJ
                IF(J.EQ.5) DUM2=DUM1*CDJ
                IF(J.EQ.8.AND.NORM) DUM2=DUM2*SQRT3
                IF(J.EQ.11) DUM2=DUM1*CFJ
                IF(J.EQ.14.AND.NORM) DUM2=DUM2*SQRT5
                IF(J.EQ.20.AND.NORM) DUM2=DUM2*SQRT3
                IF(J.EQ.21) DUM2=DUM1*CGJ
                IF(J.EQ.24.AND.NORM) DUM2=DUM2*SQRT7
                IF(J.EQ.30.AND.NORM) DUM2=DUM2*SQRT5/SQRT3
                IF(J.EQ.33.AND.NORM) DUM2=DUM2*SQRT3
C
                IJ=IJ+1
                IF(RMETHOD.EQ.RESC) THEN
                  DEN=ONE
                  IF(LOCI+I.EQ.LOCJ+J) DEN=TWO
                ELSE
                  NN = IA(MAX0(LOCI+I,LOCJ+J))+MIN0(LOCI+I,LOCJ+J)
                  DEN = DENAB(NN)
                  DEN = DEN + DEN
                  IF(RMETHOD.EQ.ANESC) THEN
                    IF(MODUS.EQ.1) DEN=DEN*FSC2
                    IF(MODUS.EQ.2) DEN=-DEN
                  ENDIF
                ENDIF
                DIJ(IJ)=DUM2*DEN
  600         CONTINUE
C
C     -----  KINETIC ENERGY
C
      IF(ICHARM.GT.0) GO TO 700
C
      T = SQRT(AA1)
      X0 = AX
      Y0 = AY
      Z0 = AZ
      DO 660 J = 1,LJTMOD
        NJ =J
        DO 660 I = 1,LITDER
          NI = I
          CALL VINT
          XS(I,J)=XINT*T
          YS(I,J)=YINT*T
          ZS(I,J)=ZINT*T
  660 CONTINUE
      CALL DTXYZ(XT,YT,ZT,XS,YS,ZS,LITDER,LJT,AJ)
      CALL DERI(DXS,DYS,DZS,XS,YS,ZS,LIT,LJT,AI)
      CALL DERI(DXT,DYT,DZT,XT,YT,ZT,LIT,LJT,AI)
      IJ=0
      DO 680 I=MINI,MAXI
        IX=IJX(I)
        IY=IJY(I)
        IZ=IJZ(I)
        DO 670 J=MINJ,MAXJ
          JX=IJX(J)
          JY=IJY(J)
          JZ=IJZ(J)
          DUMX=DXT(IX,JX)* YS(IY,JY)* ZS(IZ,JZ)
     1        +DXS(IX,JX)* YT(IY,JY)* ZS(IZ,JZ)
     2        +DXS(IX,JX)* YS(IY,JY)* ZT(IZ,JZ)
          DUMY= XT(IX,JX)*DYS(IY,JY)* ZS(IZ,JZ)
     1        + XS(IX,JX)*DYT(IY,JY)* ZS(IZ,JZ)
     2        + XS(IX,JX)*DYS(IY,JY)* ZT(IZ,JZ)
          DUMZ= XT(IX,JX)* YS(IY,JY)*DZS(IZ,JZ)
     1        + XS(IX,JX)* YT(IY,JY)*DZS(IZ,JZ)
     2        + XS(IX,JX)* YS(IY,JY)*DZT(IZ,JZ)
          IJ=IJ+1
          IF(RMETHOD.EQ.RESC) THEN
             NN=IA(MAX0(LOCI+I,LOCJ+J))+MIN0(LOCI+I,LOCJ+J)
             DUMINT=DIJ(IJ)
             DTINT(NN,1,IAT)=DTINT(NN,1,IAT)+DUMX*DUMINT
             DTINT(NN,2,IAT)=DTINT(NN,2,IAT)+DUMY*DUMINT
             DTINT(NN,3,IAT)=DTINT(NN,3,IAT)+DUMZ*DUMINT
          ELSE
             DE(1,IAT)=DE(1,IAT)+ DUMX*DIJ(IJ)
             DE(2,IAT)=DE(2,IAT)+ DUMY*DIJ(IJ)
             DE(3,IAT)=DE(3,IAT)+ DUMZ*DIJ(IJ)
          ENDIF
  670     CONTINUE
  680   CONTINUE
        IF(RMETHOD.EQ.ANESC.AND.MODUS.NE.0) GOTO 1000
C
C     ..... NUCLEAR ATTRACTION
C
  700 CONTINUE
      AAX = AA*AX
      AAY = AA*AY
      AAZ = AA*AZ
      IF(ICHARM.EQ.0) THEN
         MAXIC=NAT
         IF(ISEPS) MAXIC=NAT+NPS
      ELSE
         MAXIC=1
      END IF
C
      DO 790 IC = 1,MAXIC
         IF(ICHARM.EQ.0) THEN
            IF(IC.LE.NAT) THEN
               ZNUC = -ZAN(IC)
               CX = C(1,IC)
               CY = C(2,IC)
               CZ = C(3,IC)
            ELSE
               IAZ=IATSP(IC-NAT)
               ZNUC = -COSZAN(IC-NAT)
               CX = CORZAN(1,IC-NAT)
               CY = CORZAN(2,IC-NAT)
               CZ = CORZAN(3,IC-NAT)
            END IF
         ELSE
           IF(.NOT.ESD) THEN
            ZNUC = -QCHM(ICHARM)
            CX = XCHM(ICHARM)
            CY = YCHM(ICHARM)
            CZ = ZCHM(ICHARM)
            ENDIF
C           FOR ESD THE VARIABLES ARE SET ABOVE IN ESDDERZ.
         END IF
         XX = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
         IF(NROOTS.LE.3) CALL RT123
         IF(NROOTS.EQ.4) CALL ROOT4
         IF(NROOTS.EQ.5) CALL ROOT5
         DO 740 K = 1,NROOTS
            UU = AA*U(K)
            WW = W(K)*ZNUC
            TT = ONE/(AA+UU)
            T = SQRT(TT)
            X0 = (AAX+UU*CX)*TT
            Y0 = (AAY+UU*CY)*TT
            Z0 = (AAZ+UU*CZ)*TT
            DO 730 J = 1,LJT
               NJ = J
               DO 720 I = 1,LITDER
                  NI = I
                  CALL VINT
                  XV(I,J,K) = XINT
                  YV(I,J,K) = YINT
                  ZV(I,J,K) = ZINT*WW
  720          CONTINUE
  730       CONTINUE
            CALL DERI(DXV(1,1,K),DYV(1,1,K),DZV(1,1,K),
     *                 XV(1,1,K), YV(1,1,K), ZV(1,1,K),LIT,LJT,AI)
  740    CONTINUE
         IJ=0
         DO 780 I=MINI,MAXI
           IX=IJX(I)
           IY=IJY(I)
           IZ=IJZ(I)
           DO 780 J=MINJ,MAXJ
             JX=IJX(J)
             JY=IJY(J)
             JZ=IJZ(J)
             DUMX=ZERO
             DUMY=ZERO
             DUMZ=ZERO
             DO 770 K=1,NROOTS
               DUMX=DUMX+DXV(IX,JX,K)* YV(IY,JY,K)* ZV(IZ,JZ,K)
               DUMY=DUMY+ XV(IX,JX,K)*DYV(IY,JY,K)* ZV(IZ,JZ,K)
               DUMZ=DUMZ+ XV(IX,JX,K)* YV(IY,JY,K)*DZV(IZ,JZ,K)
  770        CONTINUE
             IJ=IJ+1
             IF((IC.GT.NAT).AND.(IAT.EQ.IAZ)) GOTO 780
             DUMINT=DIJ(IJ)*AA1*PI212
             IF(RMETHOD.EQ.RESC) THEN
                NN=IA(MAX0(LOCI+I,LOCJ+J))+MIN0(LOCI+I,LOCJ+J)
                DVINT(NN,1,IAT)=DVINT(NN,1,IAT)+DUMX*DUMINT
                DVINT(NN,2,IAT)=DVINT(NN,2,IAT)+DUMY*DUMINT
                DVINT(NN,3,IAT)=DVINT(NN,3,IAT)+DUMZ*DUMINT
C               IE. NO CHARMS NOW WITH RESC
             ELSE
                IF(ESD.AND.ICHARM.GT.0) THEN
C                  SAVE HELFEY GRADIENT TERMS TO 2ND BODY ARRAYS
                   IAZ=IXFTCH(X(LIAGLOB),IAT)
                   IOFF=LFMODE+3*NATFMO+(IAZ-1)*3
                   X(IOFF  )=X(IOFF  )+DUMX*DUMINT
                   X(IOFF+1)=X(IOFF+1)+DUMY*DUMINT
                   X(IOFF+2)=X(IOFF+2)+DUMZ*DUMINT
C                  WRITE(6,*) 'WWWADD1',IAZ,IAT
                ELSE
                   DE(1,IAT)=DE(1,IAT)+DUMX*DUMINT
                   DE(2,IAT)=DE(2,IAT)+DUMY*DUMINT
                   DE(3,IAT)=DE(3,IAT)+DUMZ*DUMINT
                ENDIF
                IF(ISEPS.AND.(IC.GT.NAT)) THEN
                   DE(1,IAZ)=DE(1,IAZ)-DUMX*(DIJ(IJ)*AA1*PI212)
                   DE(2,IAZ)=DE(2,IAZ)-DUMY*(DIJ(IJ)*AA1*PI212)
                   DE(3,IAZ)=DE(3,IAZ)-DUMZ*(DIJ(IJ)*AA1*PI212)
                END IF
                IF(ICHARM.GT.0) THEN
                   IF(ESD) THEN
C                     SAVE TVDER GRADIENT TERMS TO 2ND BODY ARRAYS
                      IOFF=LFMODE+3*NATFMO+(ICHARM-1)*3
                      X(IOFF  )=X(IOFF  )-DUMX*DUMINT
                      X(IOFF+1)=X(IOFF+1)-DUMY*DUMINT
                      X(IOFF+2)=X(IOFF+2)-DUMZ*DUMINT
C                     WRITE(6,*) 'WWWADD2',ICHARM
C            CHARMM CODE BELOW NEEDS THE FORCE, NOT THE GRADIENT?
                   ELSE
                      DXELMM(ICHARM)=DXELMM(ICHARM)+DUMX*DUMINT
                      DYELMM(ICHARM)=DYELMM(ICHARM)+DUMY*DUMINT
                      DZELMM(ICHARM)=DZELMM(ICHARM)+DUMZ*DUMINT
                   ENDIF
               END IF
             ENDIF
  780    CONTINUE
  790 CONTINUE
C
 1000 CONTINUE
 1100 CONTINUE
C
C     ----- END OF PRIMITIVE LOOPS -----
C
      IF(DBG) THEN
         WRITE(IW,9100) II,JJ
         CALL EGOUT(DE,NAT)
      END IF
 1300 CONTINUE
 1400 CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
 2000 CONTINUE
      IF(OUT) THEN
         IF(GOPARR) CALL DDI_GSUMF(1503,DE,3*NAT)
         IF(MASWRK) THEN
            WRITE(IW,9000)
            CALL EGOUT(DE,NAT)
         END IF
         IF(GOPARR) CALL DSCAL(3*NAT,ONE/NPROC,DE,1)
C              WILL ONLY PRINT PARTIAL CONTRIBUTIONS IN PARALLEL.
         IF(ESD) THEN
            IF(MASWRK) WRITE(IW,9088)
            CALL EGOUT(X(LFMODE+3*NATFMO),NATFMO)
         END IF
      END IF
C
      IF(GOPARR .AND. NXT) CALL DDI_DLBRESET
C
      IF(RMETHOD.EQ.RESC) THEN
        IF(GOPARR) CALL DDI_GSUMF(1503,DTINT,L2*3*NAT)
        IF(GOPARR) CALL DDI_GSUMF(1504,DVINT,L2*3*NAT)
        SVDSKW=DSKWRK
        DSKWRK=.TRUE.
         LREC=3*NAT
C        STARTING RECORD FOR THE T DERIVATIVE INTEGRALS
        DO I=1,NAT
          DO J=1,3
            LREC=LREC+1
            CALL RAWRIT(NQRDAF,X(MORDA),DTINT(1,J,I),L2,LREC,0)
C
C           ADD TWO CONTRIBUTIONS, HELLMANN-FEYNMAN AND DV GOT HERE
C
            CALL RAREAD(NQRDAF,X(MORDA),DTINT(1,J,I),L2,LREC+6*NAT,0)
            CALL DAXPY(L2,ONE,DTINT(1,J,I),1,DVINT(1,J,I),1)
            CALL RAWRIT(NQRDAF,X(MORDA),DVINT(1,J,I),L2,LREC+6*NAT,0)
          ENDDO
        ENDDO
        DSKWRK=SVDSKW
      ENDIF
C
      RETURN
C
 9000 FORMAT(/10X,38(1H-)/10X,'GRADIENT INCLUDING AO DERIVATIVE TERMS'/
     *        10X,38(1H-))
 9088 FORMAT(/10X,36(1H-)/10X,'ONE-ELECTRON ESD GRADIENT',/10X,36(1H-))
 9100 FORMAT(1X,'TVDER: SHELLS II,JJ=',2I5)
      END
C*MODULE GRD1    *DECK VINT
      SUBROUTINE VINT
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION H(28),W(28),MIN(7),MAX(7)
C
      COMMON /DSTV  / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,
     *                XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      COMMON /HERMIT/ H1,H2(2),H3(3),H4(4),H5(5),H6(6),H7(7),
     *                H8(8),H9(9),H10(10)
      COMMON /WERMIT/ W1,W2(2),W3(3),W4(4),W5(5),W6(6),W7(7),
     *                W8(8),W9(9),W10(10)
C
      PARAMETER (ZERO=0.0D+00)
C
      EQUIVALENCE (H(1),H1),(W(1),W1)
C
      DATA MIN /1,2,4,7,11,16,22/
      DATA MAX /1,3,6,10,15,21,28/
C
C     ----- GAUSS-HERMITE QUADRATURE USING MINIMUM POINT FORMULA -----
C
      XINT = ZERO
      YINT = ZERO
      ZINT = ZERO
      NPTS = (NI+NJ-2)/2+1
      IMIN = MIN(NPTS)
      IMAX = MAX(NPTS)
      DO 340 I = IMIN,IMAX
         DUM = W(I)
         PX = DUM
         PY = DUM
         PZ = DUM
         DUM = H(I)*T
         PTX = DUM+X0
         PTY = DUM+Y0
         PTZ = DUM+Z0
         AX = PTX-XI
         AY = PTY-YI
         AZ = PTZ-ZI
         BX = PTX-XJ
         BY = PTY-YJ
         BZ = PTZ-ZJ
         IF(PX+AX .EQ. PX) AX = ZERO
         IF(PY+AY .EQ. PY) AY = ZERO
         IF(PZ+AZ .EQ. PZ) AZ = ZERO
         IF(PX+BX .EQ. PX) BX = ZERO
         IF(PY+BY .EQ. PY) BY = ZERO
         IF(PZ+BZ .EQ. PZ) BZ = ZERO
         GO TO (180,164, 162, 160,140,120,100),NI
C
  100    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
  120    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
  140    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
  160    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
  162    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
  164    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
  180    GO TO (320,310,300,280,260,240,220,200),NJ
C
  200    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  220    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  240    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  260    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  280    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  300    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  310    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
  320    CONTINUE
         XINT = XINT+PX
         YINT = YINT+PY
         ZINT = ZINT+PZ
  340 CONTINUE
      RETURN
      END
C*MODULE GRD1    *DECK VNNDER
      SUBROUTINE VNNDER(DRG,NATM,EF3,NFRPTS,COSDRG,NDIMCD,WORK,
     *                  EFATMP,EFBTMP,CHGTMP,DIPTMP,QADTMP,
     *                  DINDTMP,DINDDTMP,EF3FRG,ATORQN,MXTMP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL OUT,GOPARR,DSKWRK,MASWRK,EFLDL,IFEWLD,TNFOIL,EWLDPL
C
      PARAMETER (MXATM=2000, MXAO=8192, MXSH=5000, MXGTOT=20000,
     *           MXDFG=5, MXPT=2000, MXFRG=1050,
     *           MXFGPT=12000, MXIFRQ=12,
     *           MXDPPT=MXFRG*MXDFG*12)
C
      DOUBLE PRECISION MD
      CHARACTER*8 POLNAM,DPOLNAM,FRGNME,DNAME,DPSNAM
C
C        MXTMP SHOULD BE JUST 1 IF EFP IS NOT IN USE,
C        OR THE GROSSLY WASTEFUL -MXFGPT- BECAUSE JIE SONG CHOSE
C        TO COPY THE ENTIRE EFP ARRAYS INSTEAD OF JUST THE PARTS
C        BEING USED...SOMEONE WHO UNDERSTANDS EFP SHOULD FIX THAT.
      DIMENSION DRG(NATM,NATM),
     *          EF3(3,*),COSDRG(NDIMCD,NDIMCD),WORK(10,*),
     *          EFATMP(MXTMP),EFBTMP(MXTMP),
     *          CHGTMP(2,MXTMP),DIPTMP(3,MXTMP),QADTMP(6,MXTMP),
     *          DINDTMP(3,MXTMP),DINDDTMP(3,MXTMP),
     *          EF3FRG(3,MXTMP),ATORQN(3,*)
      DIMENSION EF3T(3)
C
      COMMON /EFLDC / EVEC(3),EFLDL
      COMMON /EFPBAS/ EX(MXGTOT,MXDFG),CS(MXGTOT,MXDFG),
     *                CP(MXGTOT,MXDFG),CD(MXGTOT,MXDFG),
     *                CF(MXGTOT,MXDFG),CG(MXGTOT,MXDFG),
     *                PRNAME(MXFGPT),PRCORD(3,MXFGPT),EFZNUC(MXFGPT),
     *                KSTART(MXSH,MXDFG),KATOM(MXSH,MXDFG),
     *                KTYPE(MXSH,MXDFG),KNG(MXSH,MXDFG),
     *                KLOC(MXSH,MXDFG),KMIN(MXSH,MXDFG),
     *                KMAX(MXSH,MXDFG),NSHELL(MXDFG),NGAUSS(MXDFG),
     *                NATEF(MXDFG),NUMEF(MXDFG),NTPATM
      COMMON /EFPCMO/ NOCC(MXDFG),NVIR(MXDFG),NAO(MXDFG),NTAO,
     *                CANFOK(MXAO,MXDFG)
      COMMON /EFPFMO/ IEFPFMO,ISWENUC,REPNUCEFP,EFPEFP,RUNEFP
      COMMON /EFPGES/ DESFRG(6,MXFRG)
      COMMON /EFPMP2/ MP2TYP
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /EFPPAR/ EFP(3,MXFGPT),EFPOL(9,MXFGPT),
     *                EFDP(3,MXIFRQ*MXFGPT),EFDPOL(9,MXIFRQ*MXFGPT),
     *                ENO,DIND(3,MXFGPT),DINDD(3,MXFGPT),
     *                POLNAM(MXFGPT),DPOLNAM(MXIFRQ,MXFGPT),
     *                POLSCR(MXFGPT),IPLSCR,INDSCR
      COMMON /EFPOTD/ METHOD,INABIO,MOVE
      COMMON /EWALD / IFEWLD,TNFOIL,BETA,ELCNST,EPSLN,CUTOFF,CUTLST,
     *                KMAXX,NKVEC,LEVEL,EWLDPL
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FMCOM / XX(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFRG/ EFF(3,MXFRG),EFF2(3,MXFRG)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB(MXFRG),
     *                NPBF(MXFRG),NTMO
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IKREP,IEF,IP_F,NFMOPCM
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
      COMMON /POL2  / CHGHF(2,MXFGPT),DIPHF(3,MXFGPT),QADHF(6,MXFGPT),
     *                AMULTHF(MXFGPT),BMULTHF(MXFGPT),
     *                AMULT2HF(MXFGPT),BMULT2HF(MXFGPT),
     *                EFPOL2(9,MXFRG),DIND2(3,MXFGPT),DINDD2(3,MXFGPT),
     *                DCORD(3,MXPT),CDIS(3,MXFGPT),
     *                CLPD(4*MXFGPT),ZLPD(4*MXFGPT),
     *                NLPD(4*MXFGPT),KFD(MXFGPT),KLD(MXFGPT),
     *                DNAME(MXPT),DPSNAM(MXFGPT)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
C     COSMO CHANGES
C
      PARAMETER (MAXDEN=25*MXATM, LENABC=2000, NPPA=1082)
      LOGICAL ISEPS,USEPS
      COMMON /ISEPS / ISEPS, USEPS
      COMMON /SOLV  / FEPSI,RDS,DISEX2,COSURF(3,LENABC),SRAD(MXATM),
     *                QDEN(MAXDEN),AR(LENABC),
     *                NSPA,NPS,NPS2,NDEN,NPSPHER
      COMMON /SOLVI / IATSP(LENABC+1),N0(2),NP1,NP2,ISKIP
      COMMON /COSMO3/ COSZAN(NPPA),CORZAN(3,NPPA)
C
      INTEGER DDI_NP, DDI_ME
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      DATA DEBUG/8HDEBUG   /, DBUGME/8HVNNDER  /, GRD1/8HGRD1    /
      DATA MD /8HMD      /
C
C     ----- GRADIENT OF NUCLEAR REPULSION ENERGY -----
C     PLUS NUCLEAR CONTRIBUTIONS FROM VARIOUS SOLVENT MODELS
C
      CALL DDI_NPROC( DDI_NP, DDI_ME )
C
C        NOTE THAT -OUT- DOES NOT INCLUDE -MASWRK- SO ALL NODES
C        CAN CALL THE PARALLEL DEBUGGING ROUTINE (W/GLOBAL SUM).
C
      OUT = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME  .OR.
     *       EXETYP.EQ.GRD1   .OR.  NPRINT.EQ.-3)
C
C     ----- FORM DISTANCE MATRIX -----
C
      DRG(1,1) = ZERO
      DO 130 K = 2,NAT
         DRG(K,K) = ZERO
         K1 = K-1
         DO 120 L = 1,K1
            RKL = ZERO
            DO 100 I = 1,3
               RKL = RKL+(C(I,K)-C(I,L))**2
  100       CONTINUE
            DRG(K,L) = -ONE/RKL
            DRG(L,K) = SQRT(RKL)
  120    CONTINUE
  130 CONTINUE
C
C     ----- NUCLEAR REPULSION CONTRIBUTION TO GRADIENT -----
C
      DO 290 KK = 1,3
         DO 240 K = 2,NAT
            ZAK = ZAN(K)
            KM1 = K-1
            DO 230 L = 1,KM1
               ZAL = ZAN(L)
               PKL = (C(KK,K)-C(KK,L))/DRG(L,K)
               DE(KK,K) = DE(KK,K)+PKL*DRG(K,L)*ZAK*ZAL
  230       CONTINUE
  240    CONTINUE
C
         NAT1 = NAT-1
         DO 280 K = 1,NAT1
            ZAK = ZAN(K)
            KP1 = K+1
            DO 270 L = KP1,NAT
               ZAL = ZAN(L)
               PKL = (C(KK,K)-C(KK,L))/DRG(K,L)
               DE(KK,K) = DE(KK,K)+PKL*DRG(L,K)*ZAK*ZAL
  270      CONTINUE
  280    CONTINUE
  290 CONTINUE
C
C     ---- COSMO CONTRIBUTION - KIM, PUT IN BY VOLKER OCT.1998
C
      IF(ISEPS) THEN
         IF(MASWRK) WRITE(IW,*)
     *        'NOW CALCULATING SAS-NUC GRADIENT CONTRIBUTION:'
         DO 350 IZ=1,NDIMCD
            DO 355 JZ=1,NDIMCD
               COSDRG(IZ,JZ)=ZERO
  355       CONTINUE
  350    CONTINUE
C
         DO 292 K=1,NAT
            DO 294 J=1+NAT,NPS+NAT
C              IAJ=IATSP(J-NAT)
               RKJ=ZERO
               DO 296  I=1,3
                  RKJ=RKJ+(C(I,K)-CORZAN(I,J-NAT))**2
  296          CONTINUE
               COSDRG(K,J)=-1/RKJ
               COSDRG(J,K)=SQRT(RKJ)
  294       CONTINUE
  292    CONTINUE
C
      DO 297 KK=1,3
         DO  298 K=1,NAT
            ZAK=ZAN(K)
            DO 299 L=1+NAT,NPS+NAT
               IAL=IATSP(L-NAT)
               ZAL=COSZAN(L-NAT)
               PKL=(C(KK,K)-CORZAN(KK,L-NAT))/COSDRG(L,K)
               DE(KK,K)=DE(KK,K)+PKL*COSDRG(K,L)*ZAK*ZAL
               DE(KK,IAL)=DE(KK,IAL)-PKL*COSDRG(K,L)*ZAK*ZAL
  299       CONTINUE
  298    CONTINUE
  297 CONTINUE
C
      DO 360 I=1,NPS
C        IAI=IATSP(I)
         DO 362 J=1,I-1
C           IAJ=IATSP(J)
C            IF(IAI.EQ.IAJ) GOTO 362
            RIJ=ZERO
            DO 364 K=1,3
               RIJ=RIJ+(CORZAN(K,I)-CORZAN(K,J))**2
  364       CONTINUE
            COSDRG(I+NAT,J+NAT)=-1/RIJ
            COSDRG(J+NAT,I+NAT)=SQRT(RIJ)
  362    CONTINUE
  360  CONTINUE
C
      DO 366 KK=1,3
         DO 368 K=1,NPS
            IAK=IATSP(K)
            ZAK=COSZAN(K)
            DO 370 L=1,K-1
               IAL=IATSP(L)
               ZAL=COSZAN(L)
               IF(IAK.EQ.IAL) GOTO 370
               PKL=(CORZAN(KK,K)-CORZAN(KK,L))/COSDRG(L+NAT,K+NAT)
              DE(KK,IAK)=DE(KK,IAK)+PKL*COSDRG(K+NAT,L+NAT)*ZAK*ZAL
              DE(KK,IAL)=DE(KK,IAL)-PKL*COSDRG(K+NAT,L+NAT)*ZAK*ZAL
  370      CONTINUE
  368    CONTINUE
  366  CONTINUE
      ENDIF
C
C     ----- NUCLEAR CONTRIBUTIONS FROM EXTERNAL ELECTRIC FIELD -----
C
      IF (EFLDL) CALL DNUCEF(DE)
C
C     ----- NUCLEAR CONTRIBUTIONS FROM ZERNER REACTION FIELD -----
C
      IF(IZRF.EQ.1) CALL DNUCZ(DE)
C
C     ----- NUCLEAR CONTRIBUTIONS FROM EFFECTIVE FRAGMENTS -----
C
      IF(NFRPTS.GT.0) THEN
         CALL VCLR(EF3   ,1,3*NFRPTS)
         CALL VCLR(EF3FRG,1,3*NFRPTS)
      END IF
      IF(IEFC.EQ.1) CALL DNUCC(DE,EF3)
      IF(IEFD.EQ.1) CALL DNUCD(DE,EF3)
      IF(IEFQ.EQ.1) CALL DNUCQ(DE,EF3)
      IF(IEFO.EQ.1) CALL DNUCO(DE,EF3)
      IF(IEFP.EQ.1) CALL DNUCP(DE,EF3)
C---  IF(IP_F.EQ.1) CALL IEFFR(EF3)
C     - DERCCM IS CALLED IN GAMESS.SRC IF NAT>0
C       HERE IT IS CALLED FOR EFP-PCM WITHOUT QM ATOMS
      IF(IP_F.EQ.1.AND.NAT.LE.0) CALL DERCCM
C
      IF(GOPARR  .AND.  NFRG.GT.0) THEN
       DO I = 1,NFRG
         ATORQN(1,I) = ATORQ(1,I)
         ATORQN(2,I) = ATORQ(2,I)
         ATORQN(3,I) = ATORQ(3,I)
       ENDDO
       CALL VCLR(ATORQ,1,3*NFRG)
      ENDIF
C
C     ----- EFFECTIVE FRAGMENT POTENTIAL -----
C     THIS IS ALL SORTS OF FRAGMENT-FRAGMENT GRADIENT TERMS,
C     NOT NUCLEAR TERMS WHICH WERE DONE JUST ABOVE.
C
      IF(NFRG.GT.1) THEN
         IF(OUT) CALL EFP_GRD_DBG0('NUCLEAR ',EF3,ATORQN)
C
C       ----- MULTIPOLE/MULTIPOLE TERMS -----
C
C        FOR EFP/FMO, MULT/MULT TERMS ARE CALCULATED IN ANOTHER PLACE
         IF (IEFPFMO.EQ.1) THEN
            CALL DEFPCLR
            GOTO 375
         END IF
C
C            NEXT IS IMPOSSIBLE, IT IS TO KEEP THE OLD CODE FOR A WHILE.
C            THIS COMPUTATION IS NOW FOUND IN -ENUC- ROUTINE.
         IF(IW.EQ.-2341) THEN
C
         IF(IEFC.EQ.1) THEN
            CALL DCHCH
            IF(OUT) CALL EFP_GRD_DBG('DCHCH   ',DESFRG)
         END IF
         IF(IEFD.EQ.1) THEN
            CALL DCHDIP
            IF(OUT) CALL EFP_GRD_DBG('DCHDIP  ',DESFRG)
            CALL DDPDP
            IF(OUT) CALL EFP_GRD_DBG('DDPDP   ',DESFRG)
         END IF
         IF(IEFQ.EQ.1) THEN
            CALL DCHQUD(WORK)
            IF(OUT) CALL EFP_GRD_DBG('DCHQUA  ',DESFRG)
            CALL DDPQUD(WORK)
            IF(OUT) CALL EFP_GRD_DBG('DDPQUD  ',DESFRG)
            CALL DQDQD(WORK)
            IF(OUT) CALL EFP_GRD_DBG('DQDQD   ',DESFRG)
         END IF
         IF(IEFO.EQ.1) THEN
            CALL DCHOCT(WORK)
            IF(OUT) CALL EFP_GRD_DBG('DCHOCT  ',DESFRG)
         END IF
C              END OF CODE BLOCK NO LONGER IN USE.
         END IF
C
C     SINCE EFP-EFP AND EFP-FMO TERMS ARE MIXED
C     IN INDUCED DIPOLE TERMS OF EFP/FMO, YOU HAVE TO COMPUTE
C     THEM BELOW EVERY FMO CYCLE. THE REDUNDANT VALUES
C     WILL BE SUBSTRACTED LATER WITHIN MONOMER OR DIMER CALC.
C
  375    CONTINUE
C
C
C     ----- MULTIPOLE/INDUCED DIPOLE TERMS -----
C     IEFP2=0/1 IS RHF,DFT VERSUS MP2 TYPE -EFP1- POTENTIAL.
C
C     WE CALL HERE WITH EMPTY STORAGE, WHICH IS NEEDED ONLY BY EFP/FMO
C
          IFMO = 0
          CALL VALFM(LOADFM)
          LDIP1 = LOADFM + 1
          LDIP2 = LOADFM + 1
          LDFRG = LOADFM + 1
C
         IF(IEFP2.EQ.0) THEN
            IF(IEFP.EQ.1.AND.IEFC.EQ.1) THEN
               CALL DCHIND(LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DCHIND  ',DESFRG)
            END IF
            IF(IEFP.EQ.1)               THEN
               CALL DININ(LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DININ   ',DESFRG)
            END IF
            IF(IEFP.EQ.1.AND.IEFD.EQ.1) THEN
               CALL DDPIND(LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DDPIND  ',DESFRG)
            END IF
            IF(IEFP.EQ.1.AND.IEFQ.EQ.1) THEN
               CALL DQDIND(WORK,LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DQDIND  ',DESFRG)
            END IF
C
C   ---- EWALD SUMS POLARIZATION GRADIENT TERMS COME HERE
            IF (IFEWLD .AND. LEVEL .GE. 2 .AND. EWLDPL) THEN
               CALL RECIPR_POLGRAD
               CALL RECIPR_POLGRAD2
               CALL TINFOIL_POLGRAD
               CALL SELF_FIELD
               CALL EXCL_ST_POLGRAD
C--               CALL EXCL_IND_POLGRAD
C--               CALL SELF_POLGRAD
C--               CALL SELF_EXCL_POLGRAD
            END IF
         ELSE
C
C     BE CAREFUL HERE. WHEN MP2/EFP1 IS USED, THERE ARE TWO
C     LEVELS OF INDUCED DIPOLE: HF AND MP2. THEREFORE, THE
C     GRADIENT FOR DCHIND,DININ,DDPIND, AND DQDIND SHOULD
C     BE CALCULATED ON BOTH LEVELS AND ADDED TO EF3FRG AND ATORQ.
C
C     HF LEVEL FIRST:  STORE TOTAL MULTIPOLE POTENTIALS, AND LOAD
C                      THE HARTREE-FOCK LEVEL DATA ONLY.
            DO 310 M=1,MXFGPT
C     SCREEN
               EFATMP(M)=EFATRM(M)
               EFBTMP(M)=EFBTRM(M)
               EFATRM(M)=AMULTHF(M)
               EFBTRM(M)=BMULTHF(M)
C     CHARGE
               CHGTMP(1,M)=EFCHG(1,M)
               CHGTMP(2,M)=EFCHG(2,M)
               EFCHG(1,M)=CHGHF(1,M)
               EFCHG(2,M)=CHGHF(2,M)
C     DIPOLE
               DIPTMP(1,M)=EFDIP(1,M)
               DIPTMP(2,M)=EFDIP(2,M)
               DIPTMP(3,M)=EFDIP(3,M)
               EFDIP(1,M)=DIPHF(1,M)
               EFDIP(2,M)=DIPHF(2,M)
               EFDIP(3,M)=DIPHF(3,M)
C     QUADRUPLE
               QADTMP(1,M)=EFQAD(1,M)
               QADTMP(2,M)=EFQAD(2,M)
               QADTMP(3,M)=EFQAD(3,M)
               QADTMP(4,M)=EFQAD(4,M)
               QADTMP(5,M)=EFQAD(5,M)
               QADTMP(6,M)=EFQAD(6,M)
               EFQAD(1,M)=QADHF(1,M)
               EFQAD(2,M)=QADHF(2,M)
               EFQAD(3,M)=QADHF(3,M)
               EFQAD(4,M)=QADHF(4,M)
               EFQAD(5,M)=QADHF(5,M)
               EFQAD(6,M)=QADHF(6,M)
  310       CONTINUE
C
C     --USE HF LEVEL DATA: DIND,DINDD,EFCHG-EFQUD,EFATRM,EFBTRM
C
            IF(IEFP.EQ.1.AND.IEFC.EQ.1) THEN
               CALL DCHIND(LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DCHINDHF',DESFRG)
            END IF
            IF(IEFP.EQ.1)               THEN
               CALL DININ(LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DININ HF',DESFRG)
            END IF
            IF(IEFP.EQ.1.AND.IEFD.EQ.1) THEN
               CALL DDPIND(LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DDPINDHF',DESFRG)
            END IF
            IF(IEFP.EQ.1.AND.IEFQ.EQ.1) THEN
               CALL DQDIND(WORK,LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DQDINDHF',DESFRG)
            END IF
C
C    RESTORE THE TOTAL MULTIPOLE EXPANSION DATA
C
            DO M=1,MXFGPT
               EFATRM(M)=EFATMP(M)
               EFBTRM(M)=EFBTMP(M)
               EFCHG(1,M)=CHGTMP(1,M)
               EFCHG(2,M)=CHGTMP(2,M)
               EFDIP(1,M)=DIPTMP(1,M)
               EFDIP(2,M)=DIPTMP(2,M)
               EFDIP(3,M)=DIPTMP(3,M)
               EFQAD(1,M)=QADTMP(1,M)
               EFQAD(2,M)=QADTMP(2,M)
               EFQAD(3,M)=QADTMP(3,M)
               EFQAD(4,M)=QADTMP(4,M)
               EFQAD(5,M)=QADTMP(5,M)
               EFQAD(6,M)=QADTMP(6,M)
            ENDDO
C
C     MP2 LEVEL NEXT:
C     SINCE CHARGE,DIPOLE,QUADRUPLE AND EFATRM/EFBTRM
C     ARE NOW RESTORED IN COMMONS, ONLY DIND AND DINDD
C     SHOULD BE REPLACED BY DIND2 AND DINDD2 TEMPORARILY.
C
            DO M=1,MXFGPT
               DO N=1,3
                  DINDTMP(N,M) =DIND(N,M)
                  DINDDTMP(N,M)=DINDD(N,M)
                  DIND(N,M) =DIND2(N,M)
                  DINDD(N,M)=DINDD2(N,M)
               ENDDO
            ENDDO
C
            IF(IEFP.EQ.1.AND.IEFC.EQ.1) THEN
               CALL DCHIND(LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DCHINDMP',DESFRG)
            END IF
            IF(IEFP.EQ.1)               THEN
               CALL DININ(LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DININ MP',DESFRG)
            END IF
            IF(IEFP.EQ.1.AND.IEFD.EQ.1) THEN
               CALL DDPIND(LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DDPINDMP',DESFRG)
            END IF
            IF(IEFP.EQ.1.AND.IEFQ.EQ.1) THEN
               CALL DQDIND(WORK,LDIP1,LDIP2,LDFRG,IFMO)
               IF(OUT) CALL EFP_GRD_DBG('DQDINDMP',DESFRG)
            END IF
C
            DO M=1,MXFGPT
               DO N=1,3
                  DIND(N,M) =DINDTMP(N,M)
                  DINDD(N,M)=DINDDTMP(N,M)
               ENDDO
            ENDDO
         END IF
C
         IF (IEFPFMO.EQ.1) GO TO 390
C
         IF (METHOD .GE. 4) GO TO 400
C
C   ----- GRADIENT TERMS SPECIFIC TO THE -EFP1- MODEL:
C
         IF(NTMO.EQ.0) THEN
C
C                GRADIENT OF THE FITTED EXCHANGE REPULSION TERM
C
            IF(IREP.EQ.1) THEN
               MXRPTS=0
               DO IFRG=1,NFRG
                  MXRPTS = MAX(NRPTS(IFRG),MXRPTS)
               ENDDO
               MXRPTS=MXRPTS+1
               NDFRG2 = (NDFRG*NDFRG+NDFRG)/2
               LEN = MXRPTS*MXRPTS*NDFRG2
               CALL VALFM(LOADFM)
               LAREL = LOADFM + 1
               LCREL = LAREL  + LEN
               LAST  = LCREL  + LEN
               NEED = LAST - LOADFM - 1
               CALL GETFM(NEED)
               CALL DAREAD(IDAF,IODA,XX(LAREL),2*LEN,75,0)
               CALL DREPUL(XX(LAREL),XX(LCREL),MXRPTS,NDFRG2)
               IF(OUT) CALL EFP_GRD_DBG('DREPUL  ',DESFRG)
               CALL RETFM(NEED)
            END IF
C
C    THERE IS A DISPERSON TERM IN THE MP2 FLAVOR OF -EFP1- MODEL
C
            IF(IEFP2.EQ.1 .OR. MP2TYP.EQ.1) THEN
               MXDPTS=0
               DO IFRG=1,NFRG
                  MXDPTS = MAX(NDPTS(IFRG),MXDPTS)
               ENDDO
               MXDPTS=MXDPTS+1
               NDFRG2 = (NDFRG*NDFRG+NDFRG)/2
               LEN = MXDPTS*MXDPTS*NDFRG2
               CALL VALFM(LOADFM)
               LAREL = LOADFM + 1
               LDREL = LAREL  + LEN
               LCREL = LDREL  + LEN
               LXREL = LCREL  + LEN
               LBREL = LXREL  + LEN
               LEREL = LBREL  + LEN
               LAST  = LEREL  + LEN
               NEED  = LAST - LOADFM - 1
               CALL GETFM(NEED)
               CALL DAREAD(IDAF,IODA,XX(LAREL),6*LEN,81,0)
C
C     THERE'S BOTH A C6 AND C8 CONTRIBUTION TO DISPERSION
C     NOTE THAT THE -EFF2- VALUES DON'T GET PUT INTO THE GRADIENT
C     UNTIL THE CALL TO THE -EFFT- ROUTINE CALLED DURING INTEGRALS.
C
               CALL VCLR(EFF2,1,3*NFRG)
C
               CN=6.0D+00
               CALL DDISPR(XX(LAREL),XX(LDREL),XX(LCREL),
     *                     MXDPTS,NDFRG2,CN)
               IF(OUT) CALL EFP_GRD_DBG('DDISPR-6',DESFRG)
C
               CN=8.0D+00
               CALL DDISPR(XX(LXREL),XX(LBREL),XX(LEREL),
     *                     MXDPTS,NDFRG2,CN)
               IF(OUT) CALL EFP_GRD_DBG('DDISPR-8',DESFRG)
               CALL RETFM(NEED)
            END IF
         END IF
C
C        --- EFP2 EX-REP AND CHARGE TRANSFER GRADIENTS ---
C        FOR MD, THESE GRADIENTS HAVE BEEN DONE EARLIER, FOR REASONS
C        OF SPEED, SEE -GCTANA- AND -GEXREP- CALLS IN ROUTINE EFSP!
C
         IF(RUNTYP.EQ.MD) GO TO 390
         IF(NTMO.EQ.0) GO TO 390
C
C        - EXCHANGE REPULSION GRADIENTS
C
         MXMO=0
         MXBF=0
         DO I = 1, NFRG
            MXMO=MAX(MXMO,NORB(I))
            MXBF=MAX(MXBF,NPBF(I))
         END DO
         MXMO2=(MXMO*MXMO+MXMO)/2
         CALL VALFM(LOADFM)
         LPROVEC = LOADFM  + 1
         LDPVEC  = LPROVEC + MXBF*NTMO
         LFOCKMA = LDPVEC  + MXBF*NTMO*3
         LSMAT   = LFOCKMA + MXMO2*NFRG
         LTMAT   = LSMAT   + MXBF*MXBF
         LWRK    = LTMAT   + MXBF*MXBF
         LSAB    = LWRK    + MXBF
         LTAB    = LSAB    + MXMO*MXMO
         LFASQ   = LTAB    + MXMO*MXMO
         LFBSQ   = LFASQ   + MXMO*MXMO
         LDSINT  = LFBSQ   + MXMO*MXMO
         LDTINT  = LDSINT  + MXBF*MXBF*6
         LDSAB   = LDTINT  + MXBF*MXBF*6
         LDTAB   = LDSAB   + MXMO*MXMO*6
         LWRK2   = LDTAB   + MXMO*MXMO*6
         LGXRFRG = LWRK2   + MXMO*MXMO
         LXRFRG  = LGXRFRG + NFRG*NFRG*6
         LAST    = LXRFRG  + NFRG*NFRG*6*4
         NEED=LAST-LOADFM-1
         CALL GETFM(NEED)
C
         LENPV=MXBF*NTMO
         LENFM=MXMO2*NFRG
         CALL DAREAD(IDAF,IODA,XX(LPROVEC),LENPV,267,0)
         CALL DAREAD(IDAF,IODA,XX(LFOCKMA),LENFM,268,0)
         CALL DERVEC(1,XX(LDPVEC),XX(LPROVEC),MXBF,NTMO)
         CALL GEXREP(ENXR,XX(LPROVEC),XX(LFOCKMA),XX(LSMAT),
     *               XX(LTMAT),XX(LWRK),XX(LSAB),XX(LTAB),
     *               XX(LFASQ),XX(LFBSQ),XX(LDPVEC),XX(LDSINT),
     *               XX(LDTINT),XX(LDSAB),XX(LDTAB),XX(LWRK2),
     *               XX(LGXRFRG),MXBF,MXMO,XX(LXRFRG))
         CALL RETFM(NEED)
C
C        - CHARGE TRANSFER GRADIENTS
C
         IF(NTAO.EQ.0) GO TO 390
C
         MXAOEFP=0
         MXBFEFP=0
         DO I = 1, NFRG
            MXAOEFP=MAX(MXAOEFP,NAO(ISET(I)))
            MXBFEFP=MAX(MXBFEFP,NPBF(I))
         END DO
         CALL VALFM(LOADFM)
         LCANVEC = LOADFM  + 1
         LDCVEC  = LCANVEC + MXBFEFP*NTAO
         LSMAT   = LDCVEC  + MXBFEFP*NTAO*3
         LTFRG   = LSMAT   + MXBFEFP*MXBFEFP
         LTMAT   = LTFRG   + MXAOEFP*MXAOEFP*MXDFG
         LWRK    = LTMAT   + MXBFEFP*MXBFEFP
         LSAB    = LWRK    + MXBFEFP
         LTAB    = LSAB    + MXAOEFP*MXAOEFP
         LVEFP   = LTAB    + MXAOEFP*MXAOEFP
         LVABEFB = LVEFP   + MXBFEFP*MXBFEFP
         LVABEFA = LVABEFB + MXAOEFP*MXAOEFP
         LVAAEFB = LVABEFA + MXAOEFP*MXAOEFP
         LVBBEFA = LVAAEFB + MXAOEFP*MXAOEFP
         LDSINT  = LVBBEFA + MXAOEFP*MXAOEFP
         LDTINT  = LDSINT  + MXBFEFP*MXBFEFP*6
         LDSAB   = LDTINT  + MXBFEFP*MXBFEFP*6
         LDTAB   = LDSAB   + MXAOEFP*MXAOEFP*6
         LWRK2   = LDTAB   + MXAOEFP*MXAOEFP*6
         LDABEFB = LWRK2   + MXAOEFP*MXAOEFP
         LDABEFA = LDABEFB + MXAOEFP*MXAOEFP*6
         LDAAEFB = LDABEFA + MXAOEFP*MXAOEFP*6
         LDBBEFA = LDAAEFB + MXAOEFP*MXAOEFP*6
         LDVEFP  = LDBBEFA + MXAOEFP*MXAOEFP*6
         LGCTFRG = LDVEFP  + MXBFEFP*MXBFEFP*6
         LIDOIT  = LGCTFRG + NFRG*NFRG*6
         LQUAD   = LIDOIT  + NMTTPT
         LAST    = LQUAD   + 6*NMTTPT
         NEED=LAST-LOADFM-1
         CALL GETFM(NEED)
C
         LENCV=MXBFEFP*NTAO
         CALL DAREAD(IDAF,IODA,XX(LCANVEC),LENCV,269,0)
         CALL DERVEC(2,XX(LDCVEC),XX(LCANVEC),MXBFEFP,NTAO)
         CALL GCTANA(CHTNRG,XX(LCANVEC),XX(LSMAT),XX(LTFRG),
     *               XX(LTMAT),XX(LWRK),XX(LSAB),XX(LTAB),
     *               XX(LVEFP),XX(LVABEFB),XX(LVABEFA),
     *               XX(LVAAEFB),XX(LVBBEFA),XX(LDCVEC),
     *               XX(LDSINT),XX(LDTINT),XX(LDSAB),XX(LDTAB),
     *               XX(LWRK2),XX(LDABEFB),XX(LDABEFA),
     *               XX(LDAAEFB),XX(LDBBEFA),XX(LDVEFP),
     *               XX(LGCTFRG),XX(LIDOIT),XX(LQUAD),MXBFEFP,MXAOEFP)
         CALL RETFM(NEED)
C
  390    CONTINUE
      END IF
C
C         THIS ENDS ALL COMPUTATIONS OF EFP GRADIENTS,
C         EXCEPT CLEANUPS, SUMS, ETC.
C
C     ----- FRAGMENT TORQUES -----
C
  400 CONTINUE
C
      IF(GOPARR  .AND.  NFRG.GT.0) THEN
        CALL DDI_GSUMF(770,EF3FRG,3*NFRPTS)
        CALL DDI_GSUMF(776,EFF   ,3*NFRG)
        IF(IEFP2.EQ.1) CALL DDI_GSUMF(777,EFF2  ,3*NFRG)
        CALL DDI_GSUMF(774,ATORQ ,3*NFRG)
        DO I = 1,NFRG
          ATORQ(1,I) = ATORQ(1,I) + ATORQN(1,I)
          ATORQ(2,I) = ATORQ(2,I) + ATORQN(2,I)
          ATORQ(3,I) = ATORQ(3,I) + ATORQN(3,I)
        ENDDO
      ENDIF
C
      IF(IEFP2.EQ.1) THEN
        L=NMTTPT+NPTTPT+NRTTPT+NTPATM+NTMO+NDTTPT
      ELSE
        L=NMTTPT+NPTTPT+NRTTPT+NTPATM+NTMO
      END IF
C
      DO 430 I = 1,L
         DO 420 J = 1,3
            DEF(J,I) = DEF(J,I) + EF3(J,I) + EF3FRG(J,I)
C---        IF(MASWRK) WRITE(IW,*) 'I,J,DEF(J,I)',I,J,DEF(J,I)
  420    CONTINUE
  430 CONTINUE
C
      DO 440 J=1,3
         EF3T(J)=ZERO
  440 CONTINUE
      DO 480 I=1,NMTTPT+NPTTPT+NRTTPT+NTPATM+NTMO
         DO 470 J=1,3
            EF3T(J)=EF3T(J)+EF3(J,I)
  470    CONTINUE
  480 CONTINUE
C
      IF(OUT.AND.MASWRK) THEN
         WRITE(IW,9088)
         CALL EGOUT(DE,NAT)
         CALL EFOUT(EF3T,EF3T,TORQ)
      END IF
C
C     ----- DIVIDE BY NUMBER OF NODES -----
C     ALL NODES HAVE COMPUTED THE SAME QUANTITY UP TO THIS POINT,
C     SO IT MUST BE SCALED DOWN BEFORE THE EVENTUAL GLOBAL SUM.
C
      IF(GOPARR) CALL DSCAL(3*NAT,ONE/NPROC,DE,1)
      RETURN
C
 9088 FORMAT(/10X,36(1H-)/10X,'GRADIENT OF NUCLEAR REPULSION ENERGY'/
     *        10X,36(1H-))
      END
C*MODULE GRD1    *DECK GRADEN
      SUBROUTINE GRADEN(EGRAD,AU,A,U,DAB,WORK,UNCON,MODE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXATM=2000, TWO=2.0D+00)
      LOGICAL UNCON,GOPARR,DSKWRK,MASWRK,SVGPAR
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      DIMENSION DAB(*),AU(*),A(*),U(NUMU,NUM),WORK(NUM)
C
C     CONTRACT DENSITY WITH A GRADIENT MATRIX,
C     SUCH AS KINETIC DERIVATIVES.
C     IF UNCON=.T., TRANSFORM MATRIX AU (NUMU,NUMU) IN THE INTERNALLY
C     UNCONTRACTED BASIS SET INTO CONTRACTED BASIS SET,
C     USING U AS THE TRANSFORMATION MATRIX A = U-T * AU * U.
C     NOTE THAT FOR UNCON=.FALSE. A AND AU MUST BE EQUIVALENT BY CALL
C     AND MUST BE ALLOCATED DIFFERENT SPACE OTHERWISE.
C     IF MODE IS NONZERO, DIAGONAL OF A IS HALVED.
C
      L1=NUM
      IF(UNCON) THEN
C        GRADEN IS CALLED BY EACH NODE WITH DIFFERING MATRICES, THUS
C        WE CANNOT USE A PARALLELISED VERSION OF TFTRI
         SVGPAR=GOPARR
         GOPARR=.FALSE.
         CALL TFTRI(A,AU,U,WORK,L1,NUMU,NUMU)
         GOPARR=SVGPAR
      ENDIF
C
C     IF(MODE.EQ.0) THEN
C        CALL PRTRIL(A,L1)
C        CALL PRTRIL(DAB,L1)
C     ENDIF
C
C     EGRAD0=EGRAD
      IJ=0
      DO I=1,L1
        DO J=1,I
          IJ=IJ+1
          EGRAD=EGRAD+DAB(IJ)*A(IJ)*TWO
        ENDDO
C
C       HALVE THE DIAGONAL CONTRIBUTION TO COMPLY WITH THE
C       DENSITY STRUCTURE
C
        IF(MODE.NE.0) EGRAD=EGRAD-DAB(IJ)*A(IJ)
      ENDDO
      RETURN
      END
C*MODULE GRD1     *DECK EFP_GRD_DBG0
      SUBROUTINE EFP_GRD_DBG0(TAG,EF3FRG,ATORQ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      CHARACTER*8 TAG
      DIMENSION EF3FRG(3,*),ATORQ(3,*)
      LOGICAL GOPARR,MASWRK,DSKWRK
      DIMENSION BUFF(6)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C        PRINT ONLY THE FIRST POINT IN EFP GRADIENT AND TORQUE
C
      BUFF(1) = EF3FRG(1,1)
      BUFF(2) = EF3FRG(2,1)
      BUFF(3) = EF3FRG(3,1)
      BUFF(4) = ATORQ(1,1)
      BUFF(5) = ATORQ(2,1)
      BUFF(6) = ATORQ(3,1)
      CALL DDI_GSUMF(1508,BUFF,6)
      IF(MASWRK) WRITE(IW,9020) TAG,BUFF
      RETURN
 9020 FORMAT(1X,A8,': G=',3F10.7,' T=',3F10.7)
      END
C*MODULE GRD1     *DECK EFP_GRD_DBG
      SUBROUTINE EFP_GRD_DBG(TAG,GRAD_TORQ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      CHARACTER*8 TAG
      DIMENSION GRAD_TORQ(6,*)
      LOGICAL GOPARR,MASWRK,DSKWRK
      DIMENSION BUFF(6)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C        PRINT ONLY THE FIRST POINT IN EFP GRADIENT AND TORQUE
C
      BUFF(1) = GRAD_TORQ(1,2)
      BUFF(2) = GRAD_TORQ(2,2)
      BUFF(3) = GRAD_TORQ(3,2)
      BUFF(4) = GRAD_TORQ(4,2)
      BUFF(5) = GRAD_TORQ(5,2)
      BUFF(6) = GRAD_TORQ(6,2)
      CALL DDI_GSUMF(1508,BUFF,6)
      IF(MASWRK) WRITE(IW,9020) TAG,BUFF
      RETURN
 9020 FORMAT(1X,A8,': GR=',3F10.7,' TQ=',3F10.7)
      END
