C 30 APR 07 - SPW - NEW MODULE FOR NEO DERIVATIVE INTEGRALS
C                   (BASED ON GRD2B.SRC)
C
C THIS FILE CONTAINS THE CODE USED ONLY BY THE POPLE PACKAGE.
C
C*MODULE NEOG@B  *DECK JKDG80
      SUBROUTINE JKDG80_NEO(DABMAX,INEW,JNEW,KNEW,LNEW,
     *                  JTYPE,IAT,JAT,KAT,LAT,MIXED)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     ---- JKDG80 DRIVER FOR 2-E INTEGRAL DERIVATIVES FOR SP BASES ----
C     FROM GAUSSIAN-80, MICHEL DUPUIS, 6/6/85.
C     FROM HONDO8, BRETT BODE, 7/15/94.
C     REWRITTEN USING LOOPS AND INLINING FPPPP, JOSE SIERRA, 7/2001
C
      PARAMETER (MXGTOT=20000, MXSH=5000, MXATM=2000)
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      LOGICAL OUT,DBG,MIXED
      COMMON /DSHLT / RTOL,DTOL,VTOL1,VTOL2,VTOLS,OUT,DBG
CSPW  COMMON /GRAD  / FXYZ(3,MXATM)
      COMMON /GRDNEO/ FXYZ(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTAC2/ EI1,EI2,CUX
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     2                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     3                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     4                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
C
CJMS  LABELLED COMMON GSPG80 DEFINED FOR COMPUTATIONAL EFFICIENCY.
CJMS  FOR SP BASES ONLY, IT CONTAINS THE E ARRAY WHICH IS THE DAB
CJMS  ARRAY WITH INDICES IN REVERSE ORDER: E(I,J,K,L)= DAB(L,K,J,I)
CJMS  AND IS USED IN SUB JKDG80 (MODULE GRD2B). IT ORIGINATES IN:
CJMS
CJMS     1. SUBS DABCLU, DABDFT, DABGVB, DABMC, DABMP2 AND DABUMP
CJMS        (MODULE GRD2A) AND SUB DABPAU (MODULE EFPAUL) WHICH ARE
CJMS        ALL CALLED BY SUB JKDER (MODULE GRD2A)
CJMS
CJMS     2. SUB DABCLU (MODULE GRD2A) WHICH IS CALLED BY SUB EFDEN OF
CJMS        MODULE EFGRD2
CJMS
CJMS     3. SUB PAR2PDM (MODULE MP2DDI) WHICH IS CALLED BY SUB PJKDMP2
CJMS        OF MODULE MP2DDI
C
      COMMON /GSPG80/ E(0:3,0:3,0:3,0:3)
C
CJMS  LABELLED COMMON JMSG80 IS ONLY USED IN THIS MODULE GRD2B.
C
      LOGICAL IJSAME,KLSAME,IKSMJL
      COMMON /JMSG80/ PI4,PIF,P12(3,3),R12,P34(3,3),R34,FCS(4,3),E34MAX,
     2                TX21(MXG2),TH21(MXG2),TY01(MXG2),TY1Y(MXG2),
     3                TF12(MXG2),P0(3,MXG2),TE12(MXG2),CSMAB(MXG2),
     4                TX43(MXG2),TH43(MXG2),TY03(MXG2),TY3Y(MXG2),
     5                TF34(MXG2),Q0(3,MXG2),TE34(MXG2),CSMCD(MXG2),
     6                IGBEG,JGBEG,KGBEG,LGBEG,IGEND,JGEND,KGEND,LGEND,
     7                IJSAME,KLSAME,IKSMJL
C
      PARAMETER (ZER=0.0D+00, PT5=0.5D+00, ONE=1.0D+00)
CC    PARAMETER (PIDIV4=0.785398163397448309615660845819D+00)
      PARAMETER (PIDIV4=0.785398163397448D+00)
      PARAMETER (PITO52=0.349868366552497D+02)
C
CJMS  THRESHOLD VALUES TO USE 'GRDPER' INSTEAD OF THE USUAL 'GRDG80'
C
      LOGICAL   CONCPT
      DIMENSION LIMG80(3:6)
      DATA      LIMG80/ 8, 8,20,20/
C
CC    PIDIV4= ATAN(ONE)
      PI4= PIDIV4
CC    PI = PI4+PI4+PI4+PI4
CC    PITO52=(PI +PI )*PI * SQRT(PI )
      PIF= PITO52
C
C     ----- TWO ELECTRON CONTRIBUTION TO THE GRADIENT -----
C
      IJSAME= INEW.EQ.JNEW
      KLSAME= KNEW.EQ.LNEW
      IF(MIXED) THEN
         IKSMJL = .FALSE.
      ELSE
         IKSMJL=(INEW.EQ.KNEW).AND.(JNEW.EQ.LNEW)
      END IF
      IGBEG= KSTART(INEW)
      JGBEG= KSTART(JNEW)
      KGBEG= KSTART(KNEW)
      LGBEG= KSTART(LNEW)
      IGEND= IGBEG+KNG(INEW)-1
      JGEND= JGBEG+KNG(JNEW)-1
      KGEND= KGBEG+KNG(KNEW)-1
      LGEND= LGBEG+KNG(LNEW)-1
C
      MFQ= JTYPE-1
      IF(JTYPE.LE.3) MFQ= JTYPE
C
      R12= ZER
      R34= ZER
      DO M=1,3
         P12(M,1)= C(M,IAT)
         P12(M,2)= C(M,JAT)
         P12(M,3)= P12(M,2)-P12(M,1)
         R12= R12+P12(M,3)*P12(M,3)
         P34(M,1)= C(M,KAT)
         P34(M,2)= C(M,LAT)
         P34(M,3)= P34(M,2)-P34(M,1)
         R34= R34+P34(M,3)*P34(M,3)
      ENDDO
C
C     ----- PRELIMINARY SCAN OF MAGNITUDE -----
C
      LK= 0
      E34MAX=ZER
      K0= 1
      DO 125 KG=KGBEG,KGEND
         X03= EX(KG)
         CSUMC= ABS(CS(KG))+ ABS(CP(KG))
         LGND=LGEND
         IF(KLSAME) LGND=KG
         KL=K0
         DO 120 LG=LGBEG,LGND
            X04= EX(LG)
            X34= X03+X04
            X43= ONE/X34
            H43= PT5*X43
            Y03= X03*X43
            Y04= ONE-Y03
            Y34= Y03*X04
            TE34(KL)= ZER
            R34Y34= R34*Y34
            IF(R34Y34.GT.CUX) GO TO 120
            E34= X43* EXP(-R34Y34)
            IF(KLSAME.AND.(KG.NE.LG)) E34=E34+E34
            TE34(KL)= E34
            E34= E34*CSUMC*( ABS(CS(LG))+ ABS(CP(LG)))
            IF(E34MAX.LT.E34) E34MAX=E34
            CSMCD(KL)=E34*E34
            TX43(KL)= X43
            TH43(KL)= H43
            TY03(KL)= Y03
            TY3Y(KL)= Y03*Y04
            TF34(KL)= Y34+Y34
            DO 115 M=1,3
  115       Q0(M,KL)= P34(M,1)+P34(M,3)*Y04
C
            LK=LK+1
C
  120    KL=KL+1
  125 K0=K0+MXGSH
      IF(E34MAX*DABMAX.LT.VTOL1) GO TO 999
C
      JI= 0
      I0= 1
      DO 940 IG=IGBEG,IGEND
         X01= EX(IG)
         CSUMA= DABMAX*( ABS(CS(IG))+ ABS(CP(IG)))
         JGND=JGEND
         IF(IJSAME) JGND=IG
         IJ=I0
         DO 930 JG=JGBEG,JGND
            X02= EX(JG)
            X12= X01+X02
            X21= ONE/X12
            H21= PT5*X21
            Y01= X01*X21
            Y02= ONE-Y01
            Y12= Y01*X02
            TE12(IJ)= ZER
            R12Y12= R12*Y12
            IF(R12Y12.GT.CUX) GO TO 930
            E12= PIF*X21* EXP(-R12Y12)
            IF(IJSAME.AND.(IG.NE.JG)) E12=E12+E12
            CAB= E12*CSUMA*( ABS(CS(JG))+ ABS(CP(JG)))
            IF(E34MAX*CAB.LT.VTOL2) GO TO 930
            TE12(IJ)= E12
            CSMAB(IJ)=CAB*CAB
            TX21(IJ)= X21
            TH21(IJ)= H21
            TY01(IJ)= Y01
            TY1Y(IJ)= Y01*Y02
            TF12(IJ)= Y12+Y12
            DO 152 M=1,3
  152       P0(M,IJ)= P12(M,1)+P12(M,3)*Y02
C
            JI=JI+1
C
  930    IJ=IJ+1
  940 I0=I0+MXGSH
      IF(JI.EQ.0) GO TO 999
C
      DO N=1,3
         DO M=1,4
            FCS(M,N)= ZER
         ENDDO
      ENDDO
C
      CONCPT=.TRUE.
      IF(JTYPE.LT.3) THEN
C FOR JTYPE= 1 OR 2, USE THE 'GRDG80' SUBROUTINE
      ELSEIF(JTYPE.EQ.3) THEN
         IF(JI*LK.GE.LIMG80(JTYPE)) CONCPT=.FALSE.
      ELSEIF(JTYPE.EQ.4) THEN
         IF(JI*LK.GE.LIMG80(JTYPE)) CONCPT=.FALSE.
      ELSEIF(JTYPE.EQ.5) THEN
         IF(JI*LK.GE.LIMG80(JTYPE)) CONCPT=.FALSE.
      ELSEIF(JTYPE.EQ.6) THEN
         IF(JI*LK.GE.LIMG80(JTYPE)) CONCPT=.FALSE.
      ENDIF
      IF(CONCPT) THEN
         CALL GRDG80(E,JTYPE,MFQ)
      ELSE
         CALL GRDPER(E,JTYPE,MFQ)
      ENDIF
C
C SUMMATION OF THE CONTRIBUTIONS FROM THE SHELLS
C
      DO 990 M=1,3
         FCS(4,M)=-(FCS(1,M)+FCS(2,M)+FCS(3,M))
         FXYZ(M,IAT)= FXYZ(M,IAT)+FCS(1,M)
         FXYZ(M,JAT)= FXYZ(M,JAT)+FCS(2,M)
         FXYZ(M,KAT)= FXYZ(M,KAT)+FCS(3,M)
         FXYZ(M,LAT)= FXYZ(M,LAT)+FCS(4,M)
  990 CONTINUE
      IF(OUT) THEN
         WRITE(IW,9010) INEW,JNEW,KNEW,LNEW,IAT,JAT,KAT,LAT,
     2                 ((FCS(M,N),M=1,4),N=1,3)
         WRITE(IW,9020)(FXYZ(M,IAT),FXYZ(M,JAT),
     2                  FXYZ(M,KAT),FXYZ(M,LAT),M=1,3)
      ENDIF
  999 CONTINUE
C
      RETURN
 9010 FORMAT(' SHELLS AND CENTERS ',8I4/' ---- PARTIAL CONTRIBUTION'/
     2       4D15.7/4D15.7/4D15.7)
 9020 FORMAT(' ---- GRADIENT UP TO THIS POINT'/
     2       4D15.7/4D15.7/4D15.7)
      END
