C 19 SEP 05 - IA  - SYNCHRONIZE FRGINF COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL, ROOT, FRGINF, EFPPAR
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  3 JUL 03 - JMM - SUPPRESS PRINTING FOR MONTE CARLO JOBS
C 17 APR 02 - MWS - SYNCH UP EFMULT AND FRGINF COMMON
C 26 OCT 00 - MWS - INTRODUCE MXAO PARAMETER
C 16 FEB 00 - VK  - DNUC*: CHANGES DUE TO THE FROZEN QM/MM BUFFER
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C 10 SEP 96 - MWS - INCLUDE THIS CODE IN DISTRIBUTION VERSION
C 13 JUN 96 - MWS - INTRODUCE CI TYPE VARIABLE, DELETE EFPCGD
C 24 MAY 96 - WC  - EFPDER: REMOVE DAMPING TERM, AXE EFPEX IN EFPPAR
C  9 FEB 96 - MWS - DNUCP,FRGDER: DELETE SCRF CODE
C  9 JAN 96 - WC  - REMOVE EXTRANEOUS OUTPUT
C 23 NOV 94 - MWS - REMOVE ALL FTNCHECK ERRORS
C  5 JAN 95 - PND - CHANGES FOR UNSYMMETRIC POLARIZATION TENSORS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 27 MAR 92 - PND - DNUCC: ADDED ARGUMENT EF3 AND CALCULATED IT.
C 27 MAR 92 - PND - DNUCD: ADDED ARGUMENT EF3 AND CALCULATED IT.
C 27 MAR 92 - PND - DNUCQ: ADDED ARGUMENT EF3 AND CALCULATED IT.
C 27 MAR 92 - PND - DNUCO: ADDED ARGUMENT EF3 AND CALCULATED IT.
C 27 MAR 92 - PND - DNUCP: ADDED ARGUMENT EP3 AND CALCULATED IT.
C 31 OCT 91 - JHJ - EFPCGD: ADDED MCSCF OPTION.
C 28 OCT 91 - JHJ - EFPCGD: DAF FILES REPLACE FILE 25 FOR ZRF.
C 23 OCT 91 - JHJ - INTRODUCED ICALC IN COMMON /ZRFPAR/.
C 22 OCT 91 - JHJ - REPLACED OPTMIZ WITH GRAD THROUGOUT.
C 30 JUL 91 - JHJ - DNUCX(X=C,D,Q,O):ADDED DOX ARRAY.
C 23 JUL 91 - JHJ - FRGDER:INTROCUED MEMORY ALLOCATION ROUTINE.
C 22 JUL 91 - JHJ - DNUCP:DRG'S DNUCP INCLUDED.
C 17 JUL 91 - JHJ - NEW MODULE FOR FIELD CORRECTIONS TO THE GRADIENT.
C 26 FEB 91 - DRG - DIPDER FOR FIXED INDUCED FIELD GRADIENT.
C  5 JUN 90 - DRG - EFPDER FOR FIXED INDUCED MOMENT GRADIENTS.
C
C*MODULE EFGRDA  *DECK DNUCC
      SUBROUTINE DNUCC(EG3,EF3)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*8 FRGNME
C
      PARAMETER (MXATM=2000, MXAO=8192, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      LOGICAL DOMONO,DODIPO,DOQUAD,DOOCTU
C
      DIMENSION EG3(3,*),EF3(3,*)
C
      COMMON /DOMULT/ DOMONO(MXFGPT),DODIPO(MXFGPT),DOQUAD(MXFGPT),
     *                DOOCTU(MXFGPT)
      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 /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
C     EFC - NUCLEI CONTRIBUTION TO NUCLEAR GRADIENTS.
C
C
C     FROZEN BUFFER REGION EXCLUDES SOME OF THE NUCLEI... AND SOME
C     ARE INCLUDED ONLY WITH REDUCED CHARGE.
C
C     LLAT IS A COUNTER, WHILE LAT IS ATOM NUMBER
C
       IF(NBUFMO.GT.0) NATAB=NAT-LBF
       DO 305 LLAT=1,NAT
        LAT=LLAT
        ZANLAT=ZAN(LAT)
C
        IF(NBUFMO.GT.0)THEN
         IF (LLAT.GT.NATAB+NAPOL) GOTO 305
         IF (LLAT.GT.NATAB) THEN
            LAT=IAPOL(LLAT-NATAB)
            ZANLAT=POLCHG(LLAT-NATAB)
         END IF
        END IF
C
        DO 300 LEFC=1,NMTTPT
         IF(.NOT.DOMONO(LEFC))GO TO 300
         RS=(C(1,LAT)-EFC(1,LEFC))**2+(C(2,LAT)-EFC(2,LEFC))**2+
     *      (C(3,LAT)-EFC(3,LEFC))**2
         DAMP=-EFBTRM(LEFC)*EXP(-EFATRM(LEFC)*RS)
         CNSTNT=ZANLAT*EFCHG(1,LEFC)*((1+DAMP)/RS+2*EFATRM(LEFC)
     *         *DAMP)/SQRT(RS)
         CNSTNT=CNSTNT + ZANLAT*EFCHG(2,LEFC)/(RS*SQRT(RS))
         DO 295 LC=1,3
          EG3(LC,LAT)=EG3(LC,LAT)+(EFC(LC,LEFC)-C(LC,LAT))*CNSTNT
CPAUL
          EF3(LC,LEFC)=EF3(LC,LEFC)-(EFC(LC,LEFC)-C(LC,LAT))*CNSTNT
CPAUL
 295     CONTINUE
 300    CONTINUE
 305   CONTINUE
C
      RETURN
      END
C*MODULE EFGRDA  *DECK DNUCD
      SUBROUTINE DNUCD(EG3,EF3)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*8 FRGNME
C
      DIMENSION EG3(3,*),EF3(3,*),CNST3(3)
C
      PARAMETER (MXATM=2000, MXAO=8192, MXFRG=1050, MXFGPT=12000,
     *          MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      LOGICAL DOMONO,DODIPO,DOQUAD,DOOCTU
C
      COMMON /DOMULT/ DOMONO(MXFGPT),DODIPO(MXFGPT),DOQUAD(MXFGPT),
     *                DOOCTU(MXFGPT)
      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 /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
C     EFD - NUCLEI CONTRIBUTION TO NUCLEAR GRADIENTS.
C
C
C     FROZEN BUFFER REGION EXCLUDES SOME OF THE NUCLEI... AND SOME
C     ARE INCLUDED ONLY WITH REDUCED CHARGE.
C     LLAT IS A COUNTER, WHILE LAT IS ATOM NUMBER
C
       IF(NBUFMO.GT.0) NATAB=NAT-LBF
       DO 140 LLAT=1,NAT
       LAT=LLAT
       ZANLAT=ZAN(LAT)
C
       IF(NBUFMO.GT.0)THEN
         IF (LLAT.GT.NATAB+NAPOL) GOTO 140
         IF (LLAT.GT.NATAB) THEN
            LAT=IAPOL(LLAT-NATAB)
            ZANLAT=POLCHG(LLAT-NATAB)
         END IF
       END IF
C
       LEFD=0
       DO 120 IFRG=1,NFRG
       DO 120 IIEFD=1,NMPTS(IFRG)
       LEFD=LEFD+1
        IF(.NOT.DODIPO(LEFD))GO TO 120
        RS=(C(1,LAT)-EFC(1,LEFD))**2+(C(2,LAT)-EFC(2,LEFD))**2+
     *     (C(3,LAT)-EFC(3,LEFD))**2
        R3=RS*SQRT(RS)
        DAMP=EFBTRM(LEFD)*EXP(-EFATRM(LEFD)*RS)
        CNSTT1=ZANLAT*(3*(1-DAMP)/(R3*RS)-2*EFATRM(LEFD)*DAMP/R3)*
     *         (EFDIP(1,LEFD)*(C(1,LAT)-EFC(1,LEFD))+
     *          EFDIP(2,LEFD)*(C(2,LAT)-EFC(2,LEFD))+
     *          EFDIP(3,LEFD)*(C(3,LAT)-EFC(3,LEFD)))
        CNSTT2=ZANLAT*(1-DAMP)/R3
        DO 100 LC=1,3
         CNST3(LC)=CNSTT2*EFDIP(LC,LEFD)
         CNST4=CNSTT1*(EFC(LC,LEFD)-C(LC,LAT))
         CNST5=CNST3(LC)+CNST4
         EG3(LC,LAT)=EG3(LC,LAT)+CNST5
         EF3(LC,LEFD)=EF3(LC,LEFD)-CNST5
 100    CONTINUE
         DO 105 LC=1,3
         LCM=LC-1
         LC2=MOD(LCM+1,3)+1
         LC3=MOD(LCM+2,3)+1
         ATORQ(LC,IFRG)=ATORQ(LC,IFRG)+CNST3(LC2)*
     $   (C(LC3,LAT)-EFC(LC3,LEFD))
     $   -CNST3(LC3)*(C(LC2,LAT)-EFC(LC2,LEFD))
 105     CONTINUE
 120   CONTINUE
 140  CONTINUE
      RETURN
      END
C*MODULE EFGRDA  *DECK DNUCO
      SUBROUTINE DNUCO(EG3,EF3)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*8 FRGNME
C
      PARAMETER (MXATM=2000, MXAO=8192, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION EG3(3,*),EF3(3,*)
C
      LOGICAL DOMONO,DODIPO,DOQUAD,DOOCTU
C
      COMMON /DOMULT/ DOMONO(MXFGPT),DODIPO(MXFGPT),DOQUAD(MXFGPT),
     *                DOOCTU(MXFGPT)
      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 /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      DATA ONE/1.0D+00/
      DATA PT5,THREE/0.5D+00,3.0D+00/
      DATA FIVE,SIX/5.0D+00,6.0D+00/
      DATA GLIM/1.0D-14/
      DATA QLIM/1.0D-07/
C
C     EFO - NUCLEI CONTRIBUTION TO NUCLEAR GRADIENTS.
C
       I=0
       IF(NBUFMO.GT.0) NATAB=NAT-LBF
       DO 380 IFRG=1,NFRG
       DO 380 IIEFQ=1,NMPTS(IFRG)
       I=I+1
C     DO 380 I = 1,NMTTPT
      IF(.NOT.DOOCTU(I))GO TO 380
      ALFA = EFATRM(I)
      BETA = EFBTRM(I)
C
C REPLACE THIRD MOMENTS IN ARRAY SM WITH ELECTRIC OCTUPOLE TENSORS
C
      XXX2 = EFOCT(1,I)*PT5
      YYY2 = EFOCT(2,I)*PT5
      ZZZ2 = EFOCT(3,I)*PT5
      XXY2 = EFOCT(4,I)*PT5
      XXZ2 = EFOCT(5,I)*PT5
      XYY2 = EFOCT(6,I)*PT5
      YYZ2 = EFOCT(7,I)*PT5
      XZZ2 = EFOCT(8,I)*PT5
      YZZ2 = EFOCT(9,I)*PT5
      XYZ2 = EFOCT(10,I)*PT5
      DUMX = XXX2 + XYY2 + XZZ2
      DUMY = XXY2 + YYY2 + YZZ2
      DUMZ = XXZ2 + YYZ2 + ZZZ2
      XXX2 = FIVE * XXX2 - THREE * DUMX
      YYY2 = FIVE * YYY2 - THREE * DUMY
      ZZZ2 = FIVE * ZZZ2 - THREE * DUMZ
      XXY2 = FIVE * XXY2 - DUMY
      XXZ2 = FIVE * XXZ2 - DUMZ
      XYY2 = FIVE * XYY2 - DUMX
      YYZ2 = FIVE * YYZ2 - DUMZ
      XZZ2 = FIVE * XZZ2 - DUMX
      YZZ2 = FIVE * YZZ2 - DUMY
      XYZ2 = FIVE * XYZ2
C
C     FROZEN BUFFER REGION EXCLUDES SOME OF THE NUCLEI... AND SOME
C     ARE INCLUDED ONLY WITH REDUCED CHARGE.
C     JJ IS A COUNTER, WHILE J IS ATOM NUMBER
C
      DO 370 JJ = 1,NAT
       J=JJ
       ZANJ=ZAN(J)
C
       IF(NBUFMO.GT.0)THEN
         IF (JJ.GT.NATAB+NAPOL) GOTO 370
         IF (JJ.GT.NATAB) THEN
            J=IAPOL(JJ-NATAB)
            ZANJ=POLCHG(JJ-NATAB)
         END IF
       END IF
C
      X = C(1,J) - EFC(1,I)
      Y = C(2,J) - EFC(2,I)
      Z = C(3,J) - EFC(3,I)
      XX = X*X
      YY = Y*Y
      ZZ = Z*Z
      XY=X*Y
      XZ=X*Z
      YZ=Y*Z
      XXX = XX*X
      YYY = YY*Y
      ZZZ = ZZ*Z
      XXY = XX*Y
      XXZ = XX*Z
      XYY = X*YY
      YYZ = YY*Z
      XZZ = X*ZZ
      YZZ = Y*ZZ
      XYZ = X*Y*Z
      R2 = XX + YY + ZZ
      IF(R2.LT.GLIM) GO TO 370
      R7 = R2 * R2 * R2 * SQRT(R2)
      RRX = R2 * X
      RRY = R2 * Y
      RRZ = R2 * Z
      TERM = XXX2 * (FIVE*XXX - THREE*RRX)
     *     + YYY2 * (FIVE*YYY - THREE*RRY)
     *     + ZZZ2 * (FIVE*ZZZ - THREE*RRZ)
     *     + XXY2 * (FIVE*XXY - RRY)*THREE
     *     + XXZ2 * (FIVE*XXZ - RRZ)*THREE
     *     + XYY2 * (FIVE*XYY - RRX)*THREE
     *     + YYZ2 * (FIVE*YYZ - RRZ)*THREE
     *     + XZZ2 * (FIVE*XZZ - RRX)*THREE
     *     + YZZ2 * (FIVE*YZZ - RRY)*THREE
     *     + XYZ2 * (FIVE*XYZ)*SIX
      FACT = ONE
      DFACT=0
      IF(ABS(ALFA).GT.QLIM)THEN
       FACT = ONE - (BETA * EXP (-ALFA * R2))
       DFACT=2*ALFA*BETA*EXP(-ALFA*R2)
      END IF
      CNST1=ZANJ*TERM*(DFACT-7*FACT/R2)/(R7*5)
      CNSTX = ZANJ*FACT*
     * (XXX2*(15*XX)+30*XXY2*XY+30*XXZ2*XZ+
     *  XYY2*(15*YY)+XZZ2*(15*ZZ)
     *  +30*XYZ2*YZ)/(5*R7)
      CNSTXB = ZANJ*FACT*
     * (XXX2*(-6*XX)-6*YYY2*XY-6*ZZZ2*XZ-6*XXY2*XY-6*XXZ2*XZ+
     *  XYY2*(-6*XX)-6*YYZ2*XZ+XZZ2*(-6*XX)-
     *  6*YZZ2*XY)/(5*R7)
      CNSTY = ZANJ*FACT*
     * (YYY2*(15*YY)+30*YYZ2*YZ+30*XYY2*XY+
     *  YZZ2*(15*ZZ)+XXY2*(15*XX)
     *  +30*XYZ2*XZ)/(5*R7)
      CNSTYB = ZANJ*FACT*
     * (YYY2*(-6*YY)-6*ZZZ2*YZ-6*XXX2*XY-6*YYZ2*YZ-6*XYY2*XY+
     *  YZZ2*(-6*YY)-6*XZZ2*XY+XXY2*(-6*YY)-
     *  6*XXZ2*YZ)/(5*R7)
      CNSTZ = ZANJ*FACT*
     * (ZZZ2*(15*ZZ)+30*XZZ2*XZ+30*YZZ2*YZ+
     *  XXZ2*(15*XX)+YYZ2*(15*YY)
     *  +30*XYZ2*XY)/(5*R7)
      CNSTZB = ZANJ*FACT*
     * (ZZZ2*(-6*ZZ)-6*XXX2*XZ-6*YYY2*YZ-6*XZZ2*XZ-6*YZZ2*YZ+
     *  XXZ2*(-6*ZZ)-6*XXY2*YZ+YYZ2*(-6*ZZ)-
     *  6*XYY2*XZ)/(5*R7)
C
      EG3(1,J)=EG3(1,J)+(CNST1*X+CNSTX+CNSTXB)
      EG3(2,J)=EG3(2,J)+(CNST1*Y+CNSTY+CNSTYB)
      EG3(3,J)=EG3(3,J)+(CNST1*Z+CNSTZ+CNSTZB)
      EF3(1,I)=EF3(1,I)-(CNST1*X+CNSTX+CNSTXB)
      EF3(2,I)=EF3(2,I)-(CNST1*Y+CNSTY+CNSTYB)
      EF3(3,I)=EF3(3,I)-(CNST1*Z+CNSTZ+CNSTZB)
       ATORQ(1,IFRG)=ATORQ(1,IFRG)+CNSTY*Z-CNSTZ*Y
       ATORQ(2,IFRG)=ATORQ(2,IFRG)+CNSTZ*X-CNSTX*Z
       ATORQ(3,IFRG)=ATORQ(3,IFRG)+CNSTX*Y-CNSTY*X
CPAUL
 370   CONTINUE
 380   CONTINUE
      RETURN
      END
C*MODULE EFGRDA  *DECK DNUCP
      SUBROUTINE DNUCP(EG3,EF3)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*8 POLNAM,DPOLNAM
C
      DIMENSION EG3(3,*),EF3(3,*)
      DIMENSION CNST3(3)
C
      PARAMETER (MXATM=2000, MXAO=8192, MXFRG=1050,
     *         MXFGPT=12000, MXIFRQ=12,MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      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)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      DATA PT5/0.5D+00/
C
C     FIXED INDUCED DIPOLE - NUCLEI CONTRIBUTION TO NUCLEAR GRADIENTS.
C
      IF(IEFP.EQ.1)THEN
C
C     FROZEN BUFFER REGION EXCLUDES SOME OF THE NUCLEI... AND SOME
C     ARE INCLUDED ONLY WITH REDUCED CHARGE.
C     LLAT IS A COUNTER, WHILE LAT IS ATOM NUMBER
C
       IF(NBUFMO.GT.0) NATAB=NAT-LBF
       DO 140 LLAT=1,NAT
       LAT=LLAT
       ZANLAT=ZAN(LAT)
C
       IF(NBUFMO.GT.0)THEN
         IF (LLAT.GT.NATAB+NAPOL) GOTO 140
         IF (LLAT.GT.NATAB) THEN
            LAT=IAPOL(LLAT-NATAB)
            ZANLAT=POLCHG(LLAT-NATAB)
         END IF
       END IF
C
       LEFP=0
       DO 120 IFRG=1,NFRG
       DO 120 IIEFD=1,NPPTS(IFRG)
       LEFP=LEFP+1
         RS=(C(1,LAT)-EFP(1,LEFP))**2+(C(2,LAT)-EFP(2,LEFP))**2+
     *      (C(3,LAT)-EFP(3,LEFP))**2
         R3=RS*SQRT(RS)
         CNSTT1=ZANLAT*3*
     *       PT5*((DIND(1,LEFP)+DINDD(1,LEFP))*(C(1,LAT)-EFP(1,LEFP))+
     *           (DIND(2,LEFP)+DINDD(2,LEFP))*(C(2,LAT)-EFP(2,LEFP))+
     *           (DIND(3,LEFP)+DINDD(3,LEFP))*(C(3,LAT)-EFP(3,LEFP)))
     *           /(R3*RS)
C    *       ((DIND(1,LEFP))*(C(1,LAT)-EFP(1,LEFP))+
C    *           (DIND(2,LEFP))*(C(2,LAT)-EFP(2,LEFP))+
C    *           (DIND(3,LEFP))*(C(3,LAT)-EFP(3,LEFP)))
C    *           /(R3*RS)
         CNSTT2=ZANLAT/R3
C
         DO 100 LC=1,3
         CNST3(LC)=CNSTT2*PT5*(DIND(LC,LEFP)+DINDD(LC,LEFP))
C        CNST3(LC)=CNSTT2*DIND(LC,LEFP)
         CNST4=CNSTT1*(EFP(LC,LEFP)-C(LC,LAT))
         CNST5=CNST3(LC)+CNST4
         EG3(LC,LAT)=EG3(LC,LAT)+CNST5
CPAUL
          LEFPP=LEFP+NMTTPT
         EF3(LC,LEFPP)=EF3(LC,LEFPP)-CNST5
 100     CONTINUE
         DO 105 LC=1,3
         LCM=LC-1
         LC2=MOD(LCM+1,3)+1
         LC3=MOD(LCM+2,3)+1
         ATORQ(LC,IFRG)=ATORQ(LC,IFRG)+CNST3(LC2)*
     $   (C(LC3,LAT)-EFP(LC3,LEFP))
     $   -CNST3(LC3)*(C(LC2,LAT)-EFP(LC2,LEFP))
 105     CONTINUE
 120    CONTINUE
 140   CONTINUE
      END IF
C
      RETURN
      END
C*MODULE EFGRDA  *DECK DNUCQ
      SUBROUTINE DNUCQ(EG3,EF3)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*8 FRGNME
C
      LOGICAL DOMONO,DODIPO,DOQUAD,DOOCTU
C
      PARAMETER (MXATM=2000, MXAO=8192, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /DOMULT/ DOMONO(MXFGPT),DODIPO(MXFGPT),DOQUAD(MXFGPT),
     *                DOOCTU(MXFGPT)
      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 /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      DIMENSION EG3(3,*),EF3(3,*)
C
C     EFQ - NUCLEI CONTRIBUTION TO NUCLEAR GRADIENTS.
C
C     DO 55 LC=1,3
C     DO 55 IFRG=1,NFRG
C     TORQN(LC,IFRG)=ZERO
C55   CONTINUE
       IF(NBUFMO.GT.0) NATAB=NAT-LBF
       LEFQ=0
       DO 120 IFRG=1,NFRG
       DO 120 IIEFQ=1,NMPTS(IFRG)
       LEFQ=LEFQ+1
       IF(.NOT.DOQUAD(LEFQ))GO TO 120
       XX=EFQAD(1,LEFQ)
       YY=EFQAD(2,LEFQ)
       ZZ=EFQAD(3,LEFQ)
       XY=EFQAD(4,LEFQ)
       XZ=EFQAD(5,LEFQ)
       YZ=EFQAD(6,LEFQ)
       DUM=XX+YY+ZZ
       QXX=(3*XX-DUM)/2
       QYY=(3*YY-DUM)/2
       QZZ=(3*ZZ-DUM)/2
       QXY=1.5D+00*XY
       QXZ=1.5D+00*XZ
       QYZ=1.5D+00*YZ
C
C     FROZEN BUFFER REGION EXCLUDES SOME OF THE NUCLEI... AND SOME
C     ARE INCLUDED ONLY WITH REDUCED CHARGE.
C     LLAT IS A COUNTER, WHILE LAT IS ATOM NUMBER
C
       DO 100 LLAT=1,NAT
       LAT=LLAT
       ZANLAT=ZAN(LAT)
C
       IF(NBUFMO.GT.0)THEN
         IF (LLAT.GT.NATAB+NAPOL) GOTO 100
         IF (LLAT.GT.NATAB) THEN
            LAT=IAPOL(LLAT-NATAB)
            ZANLAT=POLCHG(LLAT-NATAB)
         END IF
       END IF
C
        DX=C(1,LAT)-EFC(1,LEFQ)
        DY=C(2,LAT)-EFC(2,LEFQ)
        DZ=C(3,LAT)-EFC(3,LEFQ)
        RS=DX*DX+DY*DY+DZ*DZ
        R5=RS*RS*SQRT(RS)
        R7=RS*R5
        DAMP=EFBTRM(LEFQ)*EXP(-EFATRM(LEFQ)*RS)
C       CNSTNT=ZANLAT*((2*EFATRM(LEFQ)*DAMP/R5-5*(1-DAMP)/R7)
C    &        *(QXX*(3*DX*DX-RS)+QYY*(3*DY*DY-RS)+QZZ*(3*DZ*DZ-RS)
C    &         +QXY*6*DX*DY+QXZ*6*DX*DZ+QYZ*6*DY*DZ)/3
C    &       -(1-DAMP)*(QXX*2+QYY*2+QZZ*2)/(3*R5))
        CNSTNT=ZANLAT*((2*EFATRM(LEFQ)*DAMP/R5-5*(1-DAMP)/R7)
     *        *(QXX*(3*DX*DX)+QYY*(3*DY*DY)+QZZ*(3*DZ*DZ)
     *         +QXY*6*DX*DY+QXZ*6*DX*DZ+QYZ*6*DY*DZ)/3
     *       -(1-DAMP)*(QXX*2+QYY*2+QZZ*2)/(3*R5))
        CNSTR=ZANLAT*2*(1-DAMP)/R5
        CNSTX=CNSTR*(QXX*DX+QXY*DY+QXZ*DZ)
        CNSTY=CNSTR*(QXY*DX+QYY*DY+QYZ*DZ)
        CNSTZ=CNSTR*(QXZ*DX+QYZ*DY+QZZ*DZ)
        EG3(1,LAT)=EG3(1,LAT)+CNSTNT*DX+CNSTX
        EG3(2,LAT)=EG3(2,LAT)+CNSTNT*DY+CNSTY
        EG3(3,LAT)=EG3(3,LAT)+CNSTNT*DZ+CNSTZ
        EF3(1,LEFQ)=EF3(1,LEFQ)-CNSTNT*DX-CNSTX
        EF3(2,LEFQ)=EF3(2,LEFQ)-CNSTNT*DY-CNSTY
        EF3(3,LEFQ)=EF3(3,LEFQ)-CNSTNT*DZ-CNSTZ
       ATORQ(1,IFRG)=ATORQ(1,IFRG)+CNSTY*DZ-CNSTZ*DY
       ATORQ(2,IFRG)=ATORQ(2,IFRG)+CNSTZ*DX-CNSTX*DZ
       ATORQ(3,IFRG)=ATORQ(3,IFRG)+CNSTX*DY-CNSTY*DX
 100   CONTINUE
 120  CONTINUE
      RETURN
      END
C*MODULE EFGRDA  *DECK EFPDER
      SUBROUTINE EFPDER(DM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*8 POLNAM,DPOLNAM
C
      LOGICAL ISKIP
      LOGICAL OUT,DBG,NORM,GOPARR,DSKWRK,MASWRK,NXT
C
      DIMENSION DM(*)
      DIMENSION ISKIP(35),INDX(70)
      DIMENSION IX(35),IY(35),IZ(35),JX(20),JY(20),JZ(20)
      DIMENSION IJG(210),IJX(210),IJY(210),IJZ(210)
      DIMENSION XIN(200),YIN(200),ZIN(200),G(210),DIJ(210)
      DIMENSION FIJ(210)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000,
     *           MXFRG=1050, MXFGPT=12000, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /DSTV  / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,
     *                XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      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)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,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
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)
C
      DATA CHECK/8HCHECK   /, DEBUG/8HDEBUG   /, DBUGME/8HTVDER   /
      DATA GRD1/8HGRD1    /
      DATA PT5/0.5D+00/
C
      DATA INDX / 1, 2, 3, 4, 5, 6, 7, 8, 9,10,
     1           11,12,13,14,15,16,17,18,19,20,
     2           -0,-0,-0,-0,-0,-0,-0,-0,-0,-0,
     3           -0,-0,-0,-0,-0,
     4           -0, 1, 2, 3, 1, 2, 3, 4, 5, 6,
     5            4, 5, 6, 7, 8, 9,10,11,12,13,
     6            7, 8, 9,10,11,12,13,14,15,16,
     7           17,18,19,20,21/
      DATA JX / 0, 1, 0, 0, 2, 0, 0, 1, 1, 0,
     1          3, 0, 0, 2, 2, 1, 0, 1, 0, 1/
      DATA IX / 1, 5, 1, 1, 9, 1, 1, 5, 5, 1,
     1         13, 1, 1, 9, 9, 5, 1, 5, 1, 5,
     2         17, 1, 1,13,13, 5, 1, 5, 1, 9,
     3          9, 1, 9, 5, 5/
      DATA JY / 0, 0, 1, 0, 0, 2, 0, 1, 0, 1,
     1          0, 3, 0, 1, 0, 2, 2, 0, 1, 1/
      DATA IY / 1, 1, 5, 1, 1, 9, 1, 5, 1, 5,
     1          1,13, 1, 5, 1, 9, 9, 1, 5, 5,
     2          1,17, 1, 5, 1,13,13, 1, 5, 9,
     3          1, 9, 5, 9, 5/
      DATA JZ / 0, 0, 0, 1, 0, 0, 2, 0, 1, 1,
     1          0, 0, 3, 0, 1, 0, 1, 2, 2, 1/
      DATA IZ / 1, 1, 1, 5, 1, 1, 9, 1, 5, 5,
     1          1, 1,13, 1, 5, 1, 5, 9, 9, 5,
     2          1, 1,17, 1, 5, 1, 5,13,13, 1,
     3          9, 9, 5, 5, 9/
C
C     GRADIENTS FROM FIXED INDUCED DIPOLES FOR OPTIMIZATIONS.
C
C     INITIALIZATION FOR PARALLEL
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT = -1
      LCNT = -1
      DBG = .FALSE.
      OUT = EXETYP.EQ.GRD1
      IF (EXETYP.EQ.DEBUG  .OR. EXETYP.EQ.DBUGME  .OR.
     *    NPRINT.EQ.-3) THEN
         OUT = .TRUE.
         DBG = .TRUE.
      END IF
      IF (DBG .AND. MASWRK) WRITE (IW,9008)
      TOL = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
C
      IF (EXETYP .EQ. CHECK) GO TO 1120
C
C     ----- I SHELL
C
      DO 1060 II = 1,NSHELL
C
C           GO PARALLEL!
C
      IF(NXT .AND. GOPARR) THEN
         LCNT = LCNT + 1
         IF(LCNT.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
         IF(NEXT.NE.LCNT) GO TO 1060
      END IF
      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)+1
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II)-MINI
      DO 100 I=1,35
 100  ISKIP(I)=.TRUE.
      DO 260 I=MINI,MAXI
      GO TO (110,140,260,260,180,260,260,260,260,260,
     1       220,260,260,260,260,260,260,260,260,260),I
 110  DO 120 K=2,4
 120  ISKIP(K)=.FALSE.
      GO TO 260
 140  DO 160 K=5,10
 160  ISKIP(K)=.FALSE.
      ISKIP(1)=.FALSE.
      GO TO 260
 180  DO 200 K=2,4
 200  ISKIP(K)=.FALSE.
      DO 210 K=11,20
 210  ISKIP(K)=.FALSE.
      GO TO 260
 220  DO 240 K=5,10
 240  ISKIP(K)=.FALSE.
      DO 250 K=21,35
 250  ISKIP(K)=.FALSE.
 260  CONTINUE
C
C     ----- J SHELL
C
      DO 1040 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 1040
        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
      NROOTS = (LIT+LJT+1-2)/2 + 1
      RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
C
C     ----- PREPARE INDICES FOR PAIRS OF (I,J) FUNCTIONS
C
      N0 = 0
      IF (LIT .GE. 4) N0 = 35
      IJ = 0
      DO 340 I = 1,35
      IF (ISKIP(I)) GO TO 340
      IN = INDX(I+N0)
      NX = IX(I)
      NY = IY(I)
      NZ = IZ(I)
      DO 320 J = MINJ,MAXJ
      IJ = IJ+1
      IJX(IJ) = NX+JX(J)
      IJY(IJ) = NY+JY(J)
      IJZ(IJ) = NZ+JZ(J)
      IJG(IJ) = IN+21*(J-MINJ)
  320 CONTINUE
  340 CONTINUE
      DO 360 I = 1,IJ
      N = IJG(I)
  360 G(N) = ZERO
C
C     ----- I PRIMITIVE
C
      DO 840 IG = I1,I2
      AI = EX(IG)
      ARRI = AI*RR
      AXI = AI*XI
      AYI = AI*YI
      AZI = AI*ZI
      DUM = AI+AI
      CSI = CP(IG)
      CPI = CS(IG)*DUM
      IF (LIT .EQ. 4) CPI = CD(IG)
      CDI = CP(IG)*DUM
      IF (LIT .EQ. 5) CDI = CF(IG)
      CFI = CD(IG)*DUM
      CGI = CF(IG)*DUM
C
C     ----- J PRIMITIVE
C
      DO 820 JG = J1,J2
      AJ = EX(JG)
      AA = AI+AJ
      AA1 = ONE/AA
      DUM = AJ*ARRI*AA1
      IF (DUM .GT. TOL) GO TO 820
      FAC = EXP(-DUM)
      CSJ = CS(JG)
      CPJ = CP(JG)
      CDJ = CD(JG)
      CFJ = CF(JG)
      AX = (AXI+AJ*XJ)*AA1
      AY = (AYI+AJ*YJ)*AA1
      AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
      NN = 0
      DUM1 = ZERO
      DUM2 = DUM1
      DO 600 I=1,35
      IF(ISKIP(I)) GO TO 600
      GO TO (370,380,420,420,390,420,420,420,420,420,
     1       400,420,420,420,420,420,420,420,420,420,
     2       410,420,420,420,420,420,420,420,420,420,
     3       420,420,420,420,420),I
 370  DUM1=CSI*FAC
      GO TO 420
 380  DUM1=CPI*FAC
      GO TO 420
 390  DUM1=CDI*FAC
      GO TO 420
 400  DUM1=CFI*FAC
      GO TO 420
 410  DUM1 = CGI*FAC
 420  CONTINUE
      DO 580 J=MINJ,MAXJ
      GO TO (430,440,560,560,460,560,560,480,560,560,
     1       500,560,560,520,560,560,560,560,560,540),J
  430 DUM2=DUM1*CSJ
      GO TO 560
  440 DUM2=DUM1*CPJ
      GO TO 560
  460 DUM2=DUM1*CDJ
      GO TO 560
  480 IF(NORM) DUM2=DUM2*SQRT3
      GO TO 560
  500 DUM2 = DUM1 * CFJ
      GO TO 560
  520 DUM2 = DUM2 *SQRT5
      GO TO 560
  540 DUM2 = DUM2 * SQRT3
  560 NN=NN+1
  580 DIJ(NN)=DUM2
  600 CONTINUE
C
C     ..... BASIC DIPOLE TERM
C
      DUM = 2*PI212*AA1
      DO 700 I = 1,IJ
  700 FIJ(I) = DIJ(I)*DUM
      AAX = AA*AX
      AAY = AA*AY
      AAZ = AA*AZ
      DO 819 IC = 1,NPTTPT
      CX = EFP(1,IC)
      CY = EFP(2,IC)
      CZ = EFP(3,IC)
      XMU = PT5*(DIND(1,IC)+DINDD(1,IC))
      YMU = PT5*(DIND(2,IC)+DINDD(2,IC))
      ZMU = PT5*(DIND(3,IC)+DINDD(3,IC))
      XX = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
C     CALL RT123
      IF(NROOTS.LE.3)CALL RT123
      IF(NROOTS.EQ.4)CALL ROOT4
      IF(NROOTS.EQ.5)CALL ROOT5
      MM = 0
      DO 740 K = 1,NROOTS
      UU = AA*U(K)
      WW = W(K)
      WW = WW*UU
      TT = ONE/(AA+UU)
      T = SQRT(TT)
      X0 = (AAX+UU*CX)*TT
      Y0 = (AAY+UU*CY)*TT
      Z0 = (AAZ+UU*CZ)*TT
      IN = -4+MM
      DO 720 I = 1,LIT
      IN = IN+4
      NI = I
      DO 720 J = 1,LJT
      JN = IN+J
      NJ = J
      CALL VINT
      XIN(JN) = XINT
      YIN(JN) = YINT
      ZIN(JN) = ZINT*WW
      CALL PDXYZ
      XIN(JN+100) = XINT
      YIN(JN+100) = YINT
      ZIN(JN+100) = ZINT*WW
  720 CONTINUE
  740 MM = MM+20
      DO 780 I = 1,IJ
      N = IJG(I)
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
      DUMX = ZERO
      DUMY = ZERO
      DUMZ = ZERO
      MM = 0
      DO 760 K = 1,NROOTS
      DUMX = DUMX+XIN(NX+MM+100)*YIN(NY+MM    )*ZIN(NZ+MM    )
      DUMY = DUMY+XIN(NX+MM    )*YIN(NY+MM+100)*ZIN(NZ+MM    )
      DUMZ = DUMZ+XIN(NX+MM    )*YIN(NY+MM    )*ZIN(NZ+MM+100)
  760 MM = MM+20
      DUM = FIJ(I)
  780 G(N) = G(N)-DUM*(DUMX*XMU+DUMY*YMU+DUMZ*ZMU)
 819  CONTINUE
C
C     ----- END OF PRIMITIVE LOOPS -----
C
  820 CONTINUE
  840 CONTINUE
C
C     ----- FORM INTEGRALS OVER DERIVATIVES -----
C
      NN = 0
      N = 1
      DO 900 J = MINJ,MAXJ
      IF (MINI .GT. 1) GO TO 860
      NN = NN+1
      XIN(NN) = G(N+ 1)
      YIN(NN) = G(N+ 2)
      ZIN(NN) = G(N+ 3)
      IF (MAXI .EQ. 1) GO TO 900
  860 IF (MINI .GT. 2) GO TO 880
      NN = NN+1
      XIN(NN) = (G(N+ 4)-G(N ))
      YIN(NN) = G(N+ 7)
      ZIN(NN) = G(N+ 8)
      NN = NN+1
      XIN(NN) = G(N+ 7)
      YIN(NN) = (G(N+ 5)-G(N ))
      ZIN(NN) = G(N+ 9)
      NN = NN+1
      XIN(NN) = G(N+ 8)
      YIN(NN) = G(N+ 9)
      ZIN(NN) = (G(N+ 6)-G(N ))
      IF (MAXI .EQ. 4) GO TO 900
  880 CONTINUE
      IF (MINI .GT. 5) GO TO 890
      NN = NN+1
      XIN(NN) = (G(N+ 3)-G(N )-G(N ))
      YIN(NN) = G(N+ 6)
      ZIN(NN) = G(N+ 7)
      NN = NN+1
      XIN(NN) = G(N+ 8)
      YIN(NN) = (G(N+ 4)-G(N+ 1)-G(N+ 1))
      ZIN(NN) = G(N+ 9)
      NN = NN+1
      XIN(NN) = G(N+10)
      YIN(NN) = G(N+11)
      ZIN(NN) = (G(N+ 5)-G(N+ 2)-G(N+ 2))
      NN = NN+1
      DUM = ONE
      IF (NORM) DUM = SQRT3
      XIN(NN) = DUM*(G(N+ 6)-G(N+ 1))
      YIN(NN) = DUM*(G(N+ 8)-G(N ))
      ZIN(NN) = DUM* G(N+12)
      NN = NN+1
      XIN(NN) = DUM*(G(N+ 7)-G(N+ 2))
      YIN(NN) = DUM* G(N+12)
      ZIN(NN) = DUM*(G(N+10)-G(N ))
      NN = NN+1
      XIN(NN) = DUM* G(N+12)
      YIN(NN) = DUM*(G(N+ 9)-G(N+ 2))
      ZIN(NN) = DUM*(G(N+11)-G(N+ 1))
      IF(MAXI.EQ.10) GO TO 900
 890  CONTINUE
      NN=NN+1
      XIN(NN)=(G(N+ 6)-G(N   )-G(N   )-G(N   ))
      YIN(NN)= G(N+ 9)
      ZIN(NN)= G(N+10)
      NN=NN+1
      XIN(NN)= G(N+11)
      YIN(NN)=(G(N+ 7)-G(N+ 1)-G(N+ 1)-G(N+ 1))
      ZIN(NN)= G(N+12)
      NN=NN+1
      XIN(NN)= G(N+13)
      YIN(NN)= G(N+14)
      ZIN(NN)=(G(N+ 8)-G(N+ 2)-G(N+ 2)-G(N+ 2))
      NN=NN+1
      DUM=ONE
      IF(NORM) DUM=SQRT5
      XIN(NN)=DUM* (G(N+ 9)-G(N+ 3)-G(N+ 3))
      YIN(NN)=DUM* (G(N+15)-G(N  ))
      ZIN(NN)=DUM* G(N+18)
      NN=NN+1
      XIN(NN)=DUM* (G(N+10)-G(N+ 4)-G(N+ 4))
      YIN(NN)=DUM* G(N+18)
      ZIN(NN)=DUM* (G(N+16)-G(N  ))
      NN=NN+1
      XIN(NN)=DUM* (G(N+15)-G(N+ 1))
      YIN(NN)=DUM* (G(N+11)-G(N+ 3)-G(N+ 3))
      ZIN(NN)=DUM*  G(N+19)
      NN=NN+1
      XIN(NN)=DUM* G(N+19)
      YIN(NN)=DUM* (G(N+12)-G(N+ 5)-G(N+ 5))
      ZIN(NN)=DUM* (G(N+17)-G(N+ 1))
      NN=NN+1
      XIN(NN)=DUM* (G(N+16)-G(N+ 2))
      YIN(NN)=DUM*  G(N+20)
      ZIN(NN)=DUM* (G(N+13)-G(N+ 4)-G(N+ 4))
      NN=NN+1
      XIN(NN)=DUM*  G(N+20)
      YIN(NN)=DUM* (G(N+17)-G(N+ 2))
      ZIN(NN)=DUM* (G(N+14)-G(N+ 5)-G(N+ 5))
      NN=NN+1
      IF(NORM) DUM=DUM*SQRT3
      XIN(NN)=DUM* (G(N+18)-G(N+ 5))
      YIN(NN)=DUM* (G(N+19)-G(N+ 4))
      ZIN(NN)=DUM* (G(N+20)-G(N+ 3))
  900 N = N+21
C
C     ----- CALCULATE CONTRIBUTION TO GRADIENT -----
C
      N = 0
      DO 980 J = MINJ,MAXJ
      JN = LOCJ+J
      DO 980 I = MINI,MAXI
      N = N+1
      IN = LOCI+I
      IF (IN-JN) 920,940,940
  920 ID = JN
      JD = IN
      GO TO 960
  940 ID = IN
      JD = JN
  960 NN = (ID*(ID-1))/2+JD
      DUM = DM(NN)
      DUM = DUM+DUM
      DE(1,IAT) = DE(1,IAT)+DUM*XIN(N)
      DE(2,IAT) = DE(2,IAT)+DUM*YIN(N)
      DE(3,IAT) = DE(3,IAT)+DUM*ZIN(N)
      IF (DBG .AND. MASWRK) WRITE (IW,9108) II,JJ,I,J,NN,DUM,XIN(N),
     *              YIN(N),ZIN(N),IAT,DE(1,IAT),DE(2,IAT),DE(3,IAT)
  980 CONTINUE
      IF (DBG) CALL EGOUT(DE,NAT)
 1040 CONTINUE
 1060 CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
      IF (OUT) THEN
         IF (MASWRK) WRITE (IW,9008)
         CALL EGOUT(DE,NAT)
      END IF
 1120 CONTINUE
      IF (MASWRK.AND.OUT) THEN
         WRITE (IW,9088)
         CALL TIMIT(1)
      ENDIF
      IF(GOPARR .AND. NXT) CALL DDI_DLBRESET
      RETURN
C
 9008 FORMAT(/10X,45(1H-)/10X,'EFP CONTRIBUTIONS',
     *        ' TO THE GRADIENT'/10X,45(1H-))
 9088 FORMAT(/1X,'..... END OF EFP GRADIENT ......')
 9108 FORMAT(1X,'II,JJ,I,J,NN=',4I5,I6,4E20.12/
     *       1X,'IAT,DE(1,2,3)',18X,I5,3E20.12)
      END
C*MODULE EFGRDA  *DECK FRGDER
      SUBROUTINE FRGDER(DTOT,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION DTOT(L2)
C
      PARAMETER (MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12 )
C
      COMMON /EFPOTD/ METHOF,ISCCHG,MOVE
      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 /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SIMDAT/ NACC,NREJ,IGOMIN,NRPA,IBWM,NACCT,NREJT,NRPAT,
     *                NPRTGO,IDPUNC,IGOFLG
C
      PARAMETER (ONE=1.0D+00)
      DATA CHECK/8HCHECK    /
C
C     ----- GRADIENT FOR EFFECTIVE FRAGMENT RUNS -----
C
      NGREAT = MAX(NMTTPT,NPTTPT,4*NRTTPT)
      NTOTPT = NMTTPT+NPTTPT+NRTTPT
      L4 = 300*NGREAT
C
      CALL VALFM(LOADFM)
      LCH  = LOADFM+1
      LAST = LCH + L4
      NEED = LAST - LOADFM -1
      CALL GETFM(NEED)
C
      IF(EXETYP.EQ.CHECK) GO TO 100
C
C  IF PARALLEL, SCALE NUCLEAR CONTRIBUTIONS TO FRAGMENT GRADIENT.
C   ...AND TORQUE CORRECTIONS (ATORQ)
C
      IF(GOPARR) CALL DSCAL(3*NTOTPT,ONE/NPROC,DEF,1)
      IF(GOPARR) CALL DSCAL(3*NFRG,ONE/NPROC,ATORQ,1)
C
C     ----- GET EFFECTIVE FRAGMENT DERIVATIVES -----
C
      IF(IEFC.EQ.1) CALL EFCDER(DTOT)
      IF(IEFD.EQ.1) CALL EFDDER(DTOT)
      IF(IEFQ.EQ.1) CALL EFQDER(DTOT)
      IF(IEFO.EQ.1) CALL EFODER(DTOT)
      IF(IEFP.EQ.1) CALL EFPDER(DTOT)
      IF(IREP.EQ.1) CALL REPDER(DTOT)
C
C     IF THE POSITION OF FRAGMENTS IF FIXED, SAVE CPU BU SKIPPING
C     THE GRADIENT WITH RESPECT TO FRAGMENT AND TORQUE
C
      IF(MOVE.EQ.0)GOTO 200
C
C     ----- GRADIENT WITH RESPECT TO FRAGMENT -----
C
      IF(MASWRK .AND. NPRTGO.NE.2) WRITE(IW,*)
     *      'READY TO DO 1E- PART OF GRAD WRT FRAGMENT'
      IF(IEFC.EQ.1) CALL EFCDEF(DTOT,X(LCH),L4)
      IF(IEFD.EQ.1) CALL EFDDEF(DTOT,X(LCH),L4)
      IF(IEFQ.EQ.1) CALL EFQDEF(DTOT,X(LCH),L4)
      IF(IEFO.EQ.1) CALL EFODEF(DTOT,X(LCH),L4)
      IF(IEFP.EQ.1) CALL EFPDEF(DTOT,X(LCH),L4)
      IF(IREP.EQ.1) CALL REPDEF(DTOT,X(LCH),L4)
      IF(IEFD.EQ.1) CALL EFTORD(DTOT,X(LCH),L4)
      IF(IEFQ.EQ.1) CALL EFQTOR(DTOT,X(LCH),L4)
      IF(IEFO.EQ.1) CALL EFTORO(DTOT,X(LCH),L4)
      IF(IEFP.EQ.1) CALL EFTORP(DTOT,X(LCH),L4)
  200 CONTINUE
C
C        SUM CONTRIBUTIONS FROM ALL NODES FOR PARALLEL RUNS
C
      IF(GOPARR) CALL DDI_GSUMF(1505,DEF,3*NTOTPT)
      IF(GOPARR) CALL DDI_GSUMF(1506,ATORQ,3*NFRG)
C
  100 CONTINUE
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE EFGRDA  *DECK REPDER
      SUBROUTINE REPDER(DM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*8 REPNAM
C
      LOGICAL ISKIP
      LOGICAL OUT,DBG,NORM,GOPARR,DSKWRK,MASWRK,NXT
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000,
     *           MXFRG=1050, MXFGPT=12000, 
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /DSTV  / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,
     *                XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /GRAD  / DE(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,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 /REPPAR/ CREP(3,MXFGPT),CLPR(4*MXFGPT),ZLPR(4*MXFGPT),
     *                NLPR(4*MXFGPT),KFR(MXFGPT),KLR(MXFGPT),
     *                REPNAM(MXFGPT)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
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)
C
      DIMENSION DM(2)
      DIMENSION ISKIP(35),INDX(70)
      DIMENSION IX(35),IY(35),IZ(35),JX(20),JY(20),JZ(20)
      DIMENSION IJG(210),IJX(210),IJY(210),IJZ(210)
      DIMENSION XIN(200),YIN(200),ZIN(200),G(210),DIJ(210)
      DIMENSION GIJ(210)
C
      DATA CHECK/8HCHECK   /, DEBUG/8HDEBUG   /, DBUGME/8HTVDER   /
      DATA GRD1/8HGRD1    /
      DATA INDX / 1, 2, 3, 4, 5, 6, 7, 8, 9,10,
     1           11,12,13,14,15,16,17,18,19,20,
     2           -0,-0,-0,-0,-0,-0,-0,-0,-0,-0,
     3           -0,-0,-0,-0,-0,
     4           -0, 1, 2, 3, 1, 2, 3, 4, 5, 6,
     5            4, 5, 6, 7, 8, 9,10,11,12,13,
     6            7, 8, 9,10,11,12,13,14,15,16,
     7           17,18,19,20,21/
      DATA JX / 0, 1, 0, 0, 2, 0, 0, 1, 1, 0,
     1          3, 0, 0, 2, 2, 1, 0, 1, 0, 1/
      DATA IX / 1, 5, 1, 1, 9, 1, 1, 5, 5, 1,
     1         13, 1, 1, 9, 9, 5, 1, 5, 1, 5,
     2         17, 1, 1,13,13, 5, 1, 5, 1, 9,
     3          9, 1, 9, 5, 5/
      DATA JY / 0, 0, 1, 0, 0, 2, 0, 1, 0, 1,
     1          0, 3, 0, 1, 0, 2, 2, 0, 1, 1/
      DATA IY / 1, 1, 5, 1, 1, 9, 1, 5, 1, 5,
     1          1,13, 1, 5, 1, 9, 9, 1, 5, 5,
     2          1,17, 1, 5, 1,13,13, 1, 5, 9,
     3          1, 9, 5, 9, 5/
      DATA JZ / 0, 0, 0, 1, 0, 0, 2, 0, 1, 1,
     1          0, 0, 3, 0, 1, 0, 1, 2, 2, 1/
      DATA IZ / 1, 1, 1, 5, 1, 1, 9, 1, 5, 5,
     1          1, 1,13, 1, 5, 1, 5, 9, 9, 5,
     2          1, 1,17, 1, 5, 1, 5,13,13, 1,
     3          9, 9, 5, 5, 9/
C
C
C     INITIALIZATION FOR PARALLEL
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT = -1
      LCNT = -1
      DBG = .FALSE.
      OUT = EXETYP.EQ.GRD1
      IF (EXETYP.EQ.DEBUG  .OR. EXETYP.EQ.DBUGME  .OR.
     *    NPRINT.EQ.-3) THEN
         OUT = .TRUE.
         DBG = .TRUE.
      END IF
      IF (DBG .AND. MASWRK) WRITE (IW,9008)
      TOL = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
C
      IF (EXETYP .EQ. CHECK) GO TO 1120
C
C   GRADIENT CONTRIBUTION FROM REPULSIVE POTENTIALS FOR POWERS 0
C   AND -1 OF R.
C
C     ----- I SHELL
C
      DO 1060 II = 1,NSHELL
C
C           GO PARALLEL!
C
      IF(NXT .AND. GOPARR) THEN
         LCNT = LCNT + 1
         IF(LCNT.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
         IF(NEXT.NE.LCNT) GO TO 1060
      END IF
      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)+1
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II)-MINI
      DO 100 I=1,35
 100  ISKIP(I)=.TRUE.
      DO 260 I=MINI,MAXI
      GO TO (110,140,260,260,180,260,260,260,260,260,
     1       220,260,260,260,260,260,260,260,260,260),I
 110  DO 120 K=2,4
 120  ISKIP(K)=.FALSE.
      GO TO 260
 140  DO 160 K=5,10
 160  ISKIP(K)=.FALSE.
      ISKIP(1)=.FALSE.
      GO TO 260
 180  DO 200 K=2,4
 200  ISKIP(K)=.FALSE.
      DO 210 K=11,20
 210  ISKIP(K)=.FALSE.
      GO TO 260
 220  DO 240 K=5,10
 240  ISKIP(K)=.FALSE.
      DO 250 K=21,35
 250  ISKIP(K)=.FALSE.
 260  CONTINUE
C
C     ----- J SHELL
C
      DO 1040 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 1040
        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
      NROOTS = (LIT+LJT-2)/2 + 1
      RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
C
C     ----- PREPARE INDICES FOR PAIRS OF (I,J) FUNCTIONS
C
      N0 = 0
      IF (LIT .GE. 4) N0 = 35
      IJ = 0
      DO 340 I = 1,35
      IF (ISKIP(I)) GO TO 340
      IN = INDX(I+N0)
      NX = IX(I)
      NY = IY(I)
      NZ = IZ(I)
      DO 320 J = MINJ,MAXJ
      IJ = IJ+1
      IJX(IJ) = NX+JX(J)
      IJY(IJ) = NY+JY(J)
      IJZ(IJ) = NZ+JZ(J)
      IJG(IJ) = IN+21*(J-MINJ)
  320 CONTINUE
  340 CONTINUE
      DO 360 I = 1,IJ
      N = IJG(I)
  360 G(N) = ZERO
C
C     ----- I PRIMITIVE
C
      DO 840 IG = I1,I2
      AI = EX(IG)
      ARRI = AI*RR
      AXI = AI*XI
      AYI = AI*YI
      AZI = AI*ZI
      DUM = AI+AI
      CSI = CP(IG)
      CPI = CS(IG)*DUM
      IF (LIT .EQ. 4) CPI = CD(IG)
      CDI = CP(IG)*DUM
      IF (LIT .EQ. 5) CDI = CF(IG)
      CFI = CD(IG)*DUM
      CGI = CF(IG)*DUM
C
C     ----- J PRIMITIVE
C
      DO 820 JG = J1,J2
      AJ = EX(JG)
      AA = AI+AJ
      AA1 = ONE/AA
      DUM = AJ*ARRI*AA1
      IF (DUM .GT. TOL) GO TO 820
      FAC = EXP(-DUM)
      CSJ = CS(JG)
      CPJ = CP(JG)
      CDJ = CD(JG)
      CFJ = CF(JG)
      AX = (AXI+AJ*XJ)*AA1
      AY = (AYI+AJ*YJ)*AA1
      AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
      NN = 0
      DUM1 = ZERO
      DUM2 = DUM1
      DO 600 I=1,35
      IF(ISKIP(I)) GO TO 600
      GO TO (370,380,420,420,390,420,420,420,420,420,
     1       400,420,420,420,420,420,420,420,420,420,
     2       410,420,420,420,420,420,420,420,420,420,
     3       420,420,420,420,420),I
 370  DUM1=CSI*FAC
      GO TO 420
 380  DUM1=CPI*FAC
      GO TO 420
 390  DUM1=CDI*FAC
      GO TO 420
 400  DUM1=CFI*FAC
      GO TO 420
 410  DUM1 = CGI*FAC
 420  CONTINUE
      DO 580 J=MINJ,MAXJ
      GO TO (430,440,560,560,460,560,560,480,560,560,
     1       500,560,560,520,560,560,560,560,560,540),J
  430 DUM2=DUM1*CSJ
      GO TO 560
  440 DUM2=DUM1*CPJ
      GO TO 560
  460 DUM2=DUM1*CDJ
      GO TO 560
  480 IF(NORM) DUM2=DUM2*SQRT3
      GO TO 560
  500 DUM2 = DUM1 * CFJ
      GO TO 560
  520 DUM2 = DUM2 *SQRT5
      GO TO 560
  540 DUM2 = DUM2 * SQRT3
  560 NN=NN+1
  580 DIJ(NN)=DUM2
  600 CONTINUE
      AAX = AA*AX
      AAY = AA*AY
      AAZ = AA*AZ
C
      DO 780 IC = 1,NRTTPT
      CX   =CREP(1,IC)
      CY   =CREP(2,IC)
      CZ   =CREP(3,IC)
      PCSQ = (AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2
      DO 780 LTERM=KFR(IC),KLR(IC)
      ALFA = ZLPR(LTERM)
      BETA = CLPR(LTERM)
      PREI = EXP(-AA*ALFA*PCSQ/(AA+ALFA))
      IF(NLPR(LTERM).EQ.2)THEN
C      R TO THE ZERO POWER.
       TT = ONE/(AA+ALFA)
       T  =  SQRT(TT)
       X0 = (AAX+ALFA*CX)*TT
       Y0 = (AAY+ALFA*CY)*TT
       Z0 = (AAZ+ALFA*CZ)*TT
       IN = -4
       DO 606 I = 1,LIT
       IN = IN+4
       NI = I
       DO 606 J = 1,LJT
       JN = IN+J
       NJ = J
       CALL VINT
       XIN(JN   ) = XINT*T
       YIN(JN   ) = YINT*T
       ZIN(JN   ) = ZINT*T
  606  CONTINUE
       DO 666 I = 1,IJ
       N = IJG(I)
       NX    = IJX(I)
       NY    = IJY(I)
       NZ    = IJZ(I)
       G(N)=G(N)+DIJ(I)*PREI*BETA*XIN(NX)*YIN(NY)*ZIN(NZ)
 666   CONTINUE
      ELSE
C      ONE OVER R.
       DUM = PI212/(AA+ALFA)
       DO 685 I = 1,IJ
       GIJ(I) = DIJ(I)*DUM
  685  CONTINUE
       XX=AA*AA*PCSQ/(AA+ALFA)
       IF(NROOTS .LE. 3) CALL RT123
       IF(NROOTS .EQ. 4) CALL ROOT4
       IF(NROOTS.EQ.5) CALL ROOT5
       IF(NROOTS.EQ.6) CALL ROOT6
       IF(NROOTS.GE.7)THEN
        IF (MASWRK) WRITE(IW,9008)
        CALL ABRT
       END IF
       MM = 0
       DO 725 K = 1,NROOTS
       UU = (ALFA+AA)*U(K)
       WW = W(K)
       TT = ONE/(AA+UU+ALFA)
       T  =  SQRT(TT)
       X0 = (AAX+(UU+ALFA)*CX)*TT
       Y0 = (AAY+(UU+ALFA)*CY)*TT
       Z0 = (AAZ+(UU+ALFA)*CZ)*TT
       IN = -4+MM
       DO 705 I = 1,LIT
       IN = IN+4
       NI = I
       DO 705 J = 1,LJT
       JN = IN+J
       NJ = J
       CALL VINT
       XIN(JN   ) = XINT
       YIN(JN   ) = YINT
       ZIN(JN   ) = ZINT*WW
  705  CONTINUE
  725  MM = MM+16
       DO 765 I = 1,IJ
       N = IJG(I)
       NX    = IJX(I)
       NY    = IJY(I)
       NZ    = IJZ(I)
       DUM = ZERO
       MM    = 0
       DO 745 K = 1,NROOTS
       DUM= DUM+XIN(NX+MM)*YIN(NY+MM)*ZIN(NZ+MM)
  745  MM  = MM+16
       G(N)=G(N)+GIJ(I)*PREI*BETA*DUM
 765   CONTINUE
      END IF
 780  CONTINUE
C
C     ----- END OF PRIMITIVE LOOPS -----
C
  820 CONTINUE
  840 CONTINUE
C
C     ----- FORM INTEGRALS OVER DERIVATIVES -----
C
      NN = 0
      N = 1
      DO 900 J = MINJ,MAXJ
      IF (MINI .GT. 1) GO TO 860
      NN = NN+1
      XIN(NN) = G(N+ 1)
      YIN(NN) = G(N+ 2)
      ZIN(NN) = G(N+ 3)
      IF (MAXI .EQ. 1) GO TO 900
  860 IF (MINI .GT. 2) GO TO 880
      NN = NN+1
      XIN(NN) = (G(N+ 4)-G(N ))
      YIN(NN) = G(N+ 7)
      ZIN(NN) = G(N+ 8)
      NN = NN+1
      XIN(NN) = G(N+ 7)
      YIN(NN) = (G(N+ 5)-G(N ))
      ZIN(NN) = G(N+ 9)
      NN = NN+1
      XIN(NN) = G(N+ 8)
      YIN(NN) = G(N+ 9)
      ZIN(NN) = (G(N+ 6)-G(N ))
      IF (MAXI .EQ. 4) GO TO 900
  880 CONTINUE
      IF (MINI .GT. 5) GO TO 890
      NN = NN+1
      XIN(NN) = (G(N+ 3)-G(N )-G(N ))
      YIN(NN) = G(N+ 6)
      ZIN(NN) = G(N+ 7)
      NN = NN+1
      XIN(NN) = G(N+ 8)
      YIN(NN) = (G(N+ 4)-G(N+ 1)-G(N+ 1))
      ZIN(NN) = G(N+ 9)
      NN = NN+1
      XIN(NN) = G(N+10)
      YIN(NN) = G(N+11)
      ZIN(NN) = (G(N+ 5)-G(N+ 2)-G(N+ 2))
      NN = NN+1
      DUM = ONE
      IF (NORM) DUM = SQRT3
      XIN(NN) = DUM*(G(N+ 6)-G(N+ 1))
      YIN(NN) = DUM*(G(N+ 8)-G(N ))
      ZIN(NN) = DUM* G(N+12)
      NN = NN+1
      XIN(NN) = DUM*(G(N+ 7)-G(N+ 2))
      YIN(NN) = DUM* G(N+12)
      ZIN(NN) = DUM*(G(N+10)-G(N ))
      NN = NN+1
      XIN(NN) = DUM* G(N+12)
      YIN(NN) = DUM*(G(N+ 9)-G(N+ 2))
      ZIN(NN) = DUM*(G(N+11)-G(N+ 1))
      IF(MAXI.EQ.10) GO TO 900
 890  CONTINUE
      NN=NN+1
      XIN(NN)=(G(N+ 6)-G(N   )-G(N   )-G(N   ))
      YIN(NN)= G(N+ 9)
      ZIN(NN)= G(N+10)
      NN=NN+1
      XIN(NN)= G(N+11)
      YIN(NN)=(G(N+ 7)-G(N+ 1)-G(N+ 1)-G(N+ 1))
      ZIN(NN)= G(N+12)
      NN=NN+1
      XIN(NN)= G(N+13)
      YIN(NN)= G(N+14)
      ZIN(NN)=(G(N+ 8)-G(N+ 2)-G(N+ 2)-G(N+ 2))
      NN=NN+1
      DUM=ONE
      IF(NORM) DUM=SQRT5
      XIN(NN)=DUM* (G(N+ 9)-G(N+ 3)-G(N+ 3))
      YIN(NN)=DUM* (G(N+15)-G(N  ))
      ZIN(NN)=DUM* G(N+18)
      NN=NN+1
      XIN(NN)=DUM* (G(N+10)-G(N+ 4)-G(N+ 4))
      YIN(NN)=DUM* G(N+18)
      ZIN(NN)=DUM* (G(N+16)-G(N  ))
      NN=NN+1
      XIN(NN)=DUM* (G(N+15)-G(N+ 1))
      YIN(NN)=DUM* (G(N+11)-G(N+ 3)-G(N+ 3))
      ZIN(NN)=DUM*  G(N+19)
      NN=NN+1
      XIN(NN)=DUM* G(N+19)
      YIN(NN)=DUM* (G(N+12)-G(N+ 5)-G(N+ 5))
      ZIN(NN)=DUM* (G(N+17)-G(N+ 1))
      NN=NN+1
      XIN(NN)=DUM* (G(N+16)-G(N+ 2))
      YIN(NN)=DUM*  G(N+20)
      ZIN(NN)=DUM* (G(N+13)-G(N+ 4)-G(N+ 4))
      NN=NN+1
      XIN(NN)=DUM*  G(N+20)
      YIN(NN)=DUM* (G(N+17)-G(N+ 2))
      ZIN(NN)=DUM* (G(N+14)-G(N+ 5)-G(N+ 5))
      NN=NN+1
      IF(NORM) DUM=DUM*SQRT3
      XIN(NN)=DUM* (G(N+18)-G(N+ 5))
      YIN(NN)=DUM* (G(N+19)-G(N+ 4))
      ZIN(NN)=DUM* (G(N+20)-G(N+ 3))
  900 N = N+21
C
C     ----- CALCULATE CONTRIBUTION TO GRADIENT -----
C
      N = 0
      DO 980 J = MINJ,MAXJ
      JN = LOCJ+J
      DO 980 I = MINI,MAXI
      N = N+1
      IN = LOCI+I
      IF (IN-JN) 920,940,940
  920 ID = JN
      JD = IN
      GO TO 960
  940 ID = IN
      JD = JN
  960 NN = (ID*(ID-1))/2+JD
      DUM = DM(NN)
      DUM = DUM+DUM
      DE(1,IAT) = DE(1,IAT)+DUM*XIN(N)
      DE(2,IAT) = DE(2,IAT)+DUM*YIN(N)
      DE(3,IAT) = DE(3,IAT)+DUM*ZIN(N)
      IF (DBG .AND. MASWRK) WRITE (IW,9108) II,JJ,I,J,NN,DUM,XIN(N),
     *              YIN(N),ZIN(N),IAT,DE(1,IAT),DE(2,IAT),DE(3,IAT)
  980 CONTINUE
      IF (DBG) CALL EGOUT(DE,NAT)
 1040 CONTINUE
 1060 CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
      IF (OUT) THEN
         IF (MASWRK) WRITE (IW,9008)
         CALL EGOUT(DE,NAT)
      END IF
 1120 CONTINUE
      IF (MASWRK.AND.OUT) THEN
         WRITE (IW,9088)
         CALL TIMIT(1)
      ENDIF
      IF(GOPARR .AND. NXT) CALL DDI_DLBRESET
      RETURN
C
 9008 FORMAT(/10X,45(1H-)/10X,'EFR CONTRIBUTIONS',
     *        ' TO THE GRADIENT'/10X,45(1H-))
 9088 FORMAT(/1X,'..... END OF EFR GRADIENT ......')
 9108 FORMAT(1X,'II,JJ,I,J,NN=',4I5,I6,4E20.12/
     *       1X,'IAT,DE(1,2,3)',18X,I5,3E20.12)
      END
