C 11 APR 08 - MWS - RECUR: FIX CODE FOR LMAX=3 BY PROVIDING ALL CASES
C 19 SEP 05 - MWS - ADD TRUE NUCLEAR CHARGE ARRAY TO INFOA COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 26 MAR 03 - JMS - EXTENSIVE REWRITE
C 26 MAR 02 - KRG - USE ABRT CALL
C 20 FEB 01 - CGP - ECPINT: FIX FOR MOROKUMA DECOMPOSITION
C 26 OCT 00 - MWS - INTRODUCE MXAO PARAMETER
C 13 MAR 99 - MWS - ECPINT: FIX LENGTH OF INTEGER DAF RECORDS
C 31 DEC 98 - MWS - ECPA14: CHANGE ZEROING OF CKLU ARRAY
C  1 DEC 98 - BMB - COMPLETE REWRITE OF THE ECP CODE
C 12 NOV 98 - GDF - ELIMINATE USE OF ISOOUT
C 11 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C  6 AUG 93 - BMB - ECPINT,ECPDER: USE SYMMETRY, CHANGE USE OF MEMORY
C 22 JUN 93 - BMB - MODIFIED FOR F-FUNCTIONS, THANKS TRC
C  5 NOV 92 - TLW - MOVE THE NXTVAL LOOP
C 20 MAR 92 - MWS - CLEAN UP COMMONS /COD*EC/
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C  6 JAN 92 - TLW - ECDINT: MAKE PARALLEL AND CLEAN UP WRITING
C  4 JAN 92 - TLW - MAKE WRITES PARALLEL;ADD COMMON PAR
C 20 JUN 91 - MWS - FJ: SAVE DM VARIABLE
C  7 AUG 90 - TLW - ADD CF AND CG TO COMMON NSHEL
C  4 JAN 90 - STE - FJ,FJECP,FM,FSI: SAVE STATMENTS ADDED
C 16 DEC 89 - MWS - ADD /IOFILE/ TO ECPDER.
C 28 NOV 89 - KAN - IMPLEMENT ECP IN GAMESS FROM HONDO 7.
C
C  A BRIEF EXPLANATION OF SOME OF THE 'MAGIC' NUMBERS USED IN THIS
C  FILE AND IN ECPLIB.SRC.  FIRST THE NUMBERS ORIGINATE FROM THE
C  NUMBER OF FUNCTIONS IN EACH SHELL, I.E. S->1 P->3 D->6 F->10 G->15
C  SOME MAGIC NUMBERS ARE:  S+P+D+F -> 20
C                           S+P+D+F+G -> 35
C  AND FOR DERIVATIVES: <D|H|D'> = <D|H|P> + <D|H|F> = 6(3+10) = 78
C                        <F|H|F'> = <F|H|D> + <F|H|F> = 10(6+15)= 210
C  -NCOEF1- AND -NCOEF2- DEPEND ON THE HIGHEST ANGULAR MOMENTUM
C  IN THE BASIS SET, AND ARE DETERMINED IN ECPLIB.SRC
C
C*MODULE ECP     *DECK DAWERF
      DOUBLE PRECISION FUNCTION DAWERF(Y)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /DERFCM/ C(246),IFRST(40),ILAST(40),H
C
C     -----  ROUTINE EVALUATES THE HYBRID DAWSON-ERROR FUNCTION  -----
C
      PARAMETER (ZP5=0.5D+00, ONE=1.0D+00, TEN=10.0D+00)
C
      X= ABS(Y)
      IF(X.GE.TEN) GO TO 115
C
      XN = X/H
      NX = INT(XN)+1
      I1 = IFRST(NX)
      IL = ILAST(NX)
      T= C(IL)
      DO 110 K=IL-1,I1,-1
  110 T= C(K)+T*X
      DAWERF= T
      RETURN
C
C        LARGE X EXPANSION
C        THE ORIGINAL PROGRAM HAD AN IL=4 EXPANSION, SHOWN BELOW.
C        THE RELATIVE ACCURACY IS IMPROVED FROM 2.E-10 TO 7.E-15
C
  115 CONTINUE
      X2I = ONE/(X*X)
      IL = 8
      TXT=(IL +ZP5)*X2I
      SUM= TXT
      DO 120 K=IL,1,-1
         TXT= TXT-X2I
  120 SUM= TXT*(ONE+SUM)
      DAWERF= X*TXT*(ONE+SUM)
      RETURN
C
C---    X2I = ONE/(X*X)
C---    TXT = X2I*ZP5
C---    TX3 = X2I+TXT
C---    TX5 = X2I+TX3
C---    TX7 = X2I+TX5
C---    TX9 = X2I+TX7
C---    DAWERF= X*TXT*(ONE+TXT*(ONE+TX3*(ONE+TX5*(ONE+TX7*(ONE+TX9)))))
C
      END
C*MODULE ECP     *DECK DAWF
      DOUBLE PRECISION FUNCTION DAWF(Y)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /DAWFCM/ C(249),IFRST(40),ILAST(40),H
C
C        -----  ROUTINE EVALUATES THE DAWSON FUNCTION -----
C        JOSE SIERRA HAS TESTED THIS ROUTINE AGAINST THE NAG LIBRARY,
C        AND FINDS RELATIVE DISCREPANCIES LIE BELOW 1.E-10.  THE
C        INTERVAL H IS SMALLER HERE, SO JOSE THINKS THIS ROUTINE
C        MAY BE THE MORE CORRECT ONE.
C
      PARAMETER (ZER=0.0D+00, ZP5=0.5D+00, ONE=1.0D+00, TEN=10.0D+00)
C
      X= Y
      S= ONE
      IF(X.GE.ZER) GO TO 105
      X=-Y
      S=-ONE
  105 IF(X.GE.TEN) GO TO 115
      XN = X/H
      NX = INT(XN)+1
      I1 = IFRST(NX)
      IL = ILAST(NX)
      T= C(IL)
      DO 110 K=IL-1,I1,-1
  110 T= C(K)+T*X
      DAWF  = T*S
      RETURN
C
C        LARGE X EXPANSION
C        THE ORIGINAL PROGRAM HAD AN IL=3 EXPANSION, SHOWN BELOW.
C        THE RELATIVE ACCURACY IS IMPROVED FROM 3.E-9 TO 7.E-15
C
  115 CONTINUE
      X2I = ONE/(X*X)
      IL = 8
      TXT=(IL +ZP5)*X2I
      SUM= TXT
      DO 120 K=IL,1,-1
         TXT= TXT-X2I
  120 SUM= TXT*(ONE+SUM)
      DAWF = Y*TXT*(ONE+SUM)
      RETURN
C
C---    X2I = ONE/(X*X)
C---    TXT = X2I*ZP5
C---    TX3 = X2I+TXT
C---    TX5 = X2I+TX3
C---    TX7 = X2I+TX5
C---    DAWF  = Y*TXT*(ONE+TXT*(ONE+TX3*(ONE+TX5*(ONE+TX7))))
C
      END
C*MODULE ECP     *DECK ECPAA1
      SUBROUTINE ECPAA1(NPNP,FPQR,COEFI,COEFJ,DCOEF4,G)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C ROUTINE ECPAA1 GATHERS ROUTINES ECPA11 AND ECPA21
C
      DIMENSION FPQR(25,25,25),COEFI(*),COEFJ(*),DCOEF4(*),G(*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
C ... ECP PARAMETERS ... GENERATED IN ECPPAR
C
      COMMON /ECP1  / X01,CAX,CAY,CAZ,CA,XCA,YCA,ZCA,
     *                X02,BAX,BAY,BAZ,BA,XBA,YBA,ZBA,
     *                PHASE,DAX,DAY,DAZ,DA,XDA,YDA,ZDA,XINT,KCNTR
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /ECPDIM/ NCOEF1,NCOEF2,J1LEN,J2LEN,LLIM,NLIM,NTLIM,J4LEN
C IAMIN/MAX ARE THE ANGULAR   FUNCTION MIN AND MAX FOR I
C IPMIN/MAX ARE THE PRIMITIVE FUNCTION MIN AND MAX FOR I
C KF1/KL1 ARE THE FIRST AND LAST PRIMITIVES FORMING THE LMAX ECP
C     AND THUS THE TYPE 1 INTEGRAL
      LOGICAL IANDJ,NORM,NORMI,NORMJ
      COMMON /ECPIDX/ Q2,IAMIN,IAMAX,JAMIN,JAMAX,IPMIN,IPMAX,JPMIN,
     *                JPMAX,KF1,KL1,LLMX,NPC,NPB,IANDJ,NORM,NORMI,NORMJ
C
      COMMON /GBASE / NFST(8),NX(84),NY(84),NZ(84)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
C
      PARAMETER (ZER=0.0D+00, ONE=1.0D+00)
      PARAMETER (SR3=1.73205080756888D+00, SR5=2.23606797749979D+00)
      PARAMETER (SR7=2.64575131106459D+00)
      PARAMETER (S35=SR3*SR5, S57=SR5*SR7, S53=S57/SR3)
C
      DIMENSION ECPFAC(35)
      DATA ECPFAC/7*ONE,3*SR3,3*ONE,6*SR5,1*S35,3*ONE,6*SR7,3*S53,3*S57/
C
C ROUTINE ECPA11
C
C ECP TYPE 1 INTEGRALS FOR THE SPECIAL CASE OF ALL
C THREE FUNCTIONS ON THE SAME CENTER <AAA> (IE. ICAB = 1)
C     G WILL CONTAIN THE GRADIENT VECTOR UPON EXIT (NOT SYMMETRIZED)
C     NPNP = N+N' GIVES THE MAXIMUM ANGULAR MOMENTUM (KAPPA) NEEDED
C NOTE: L SHELLS CAN NOT BE DONE WITH THIS ROUTINE DUE TO THE COEFF
C       OF THE PRIMS MULTIPLIED IN THE RADIAL INTEGRAL SECTION
C
C THE TOTAL N FOR THE TWO SHELLS MUST BE EVEN OR THE INTEGRAL IS
C IDENTICALLY ZERO
C
      NBC = NPNP-1
      IF(MOD(NBC,2).EQ.1) RETURN
C
C LOOP OVER THE PRIMITIVES SUMMING THE RADIAL INTEGRALS IN FPSUM
C
      FPSUM = ZER
      DO 120 IG=IPMIN,IPMAX
         X01= EX(IG)
         DO 120 JG=JPMIN,JPMAX
            X02= EX(JG)
            X12= X01+X02
            DUM= COEFI(IG-IPMIN+1)*COEFJ(JG-JPMIN+1)
            FPTEMP= ZER
            DO 110 K=KF1,KL1
  110       FPTEMP= FPTEMP+FA(NBC+NLP(K),X12+ZLP(K))*CLP(K)
  120 FPSUM = FPSUM+FPTEMP*DUM
C
C NOW RETRIEVE THE ANGULAR INTEGRAL AND MULTIPLY IT BY THE
C RADIAL INTEGRAL SUM AND STORE IN G
C
      FPSQ2 = FPSUM*Q2
      MMAX= JAMAX
      NN= 1
      DO 140 I=IAMIN,IAMAX
         DUMI= FPSQ2
         IF(NORMI) DUMI= DUMI*ECPFAC(I)
         IF(IANDJ) MMAX= I
         DO 140 J=JAMIN,MMAX
            DUMJ= DUMI
            IF(NORMJ) DUMJ= DUMJ*ECPFAC(J)
            MX= NX(I)+NX(J)+1
            MY= NY(I)+NY(J)+1
            MZ= NZ(I)+NZ(J)+1
            G(NN)= G(NN)+DUMJ*FPQR(MX,MY,MZ)
  140 NN= NN+1
C
      IF(LLMX.LT.2) RETURN
C
C ROUTINE ECPA21
C
C ECP TYPE 2 INTEGRALS FOR THE SPECIAL CASE OF ALL
C THREE FUNCTIONS ON THE SAME CENTER <AAA> (IE. ICAB = 1)
C     G WILL CONTAIN THE GRADIENT VECTOR UPON EXIT (NOT SYMMETRIZED)
C     NPNP = N+N' GIVES THE MAXIMUM ANGULAR MOMENTUM (KAPPA) NEEDED
C NOTE: L SHELLS CAN NOT BE DONE WITH THIS ROUTINE DUE TO THE COEFF
C       OF THE PRIMS MULTIPLIED IN THE RADIAL INTEGRAL SECTION
C
C LOOP OVER THE ECP POTENTIALS ON CENTER K
C
      DO 300 LL=2,LLMX
         L = LL-2
         NLM1=(L*(L+1))
         KF = KFRST(KCNTR,LL)
         KL = KLAST(KCNTR,LL)
C
C LOOP OVER THE PRIMITIVES SUMMING THE RADIAL INTEGRALS IN FPSUM
C
         FPSUM = ZER
         DO 220 IG=IPMIN,IPMAX
            X01= EX(IG)
            DO 220 JG=JPMIN,JPMAX
               X02= EX(JG)
               X12= X01+X02
               DUM= COEFI(IG-IPMIN+1)*COEFJ(JG-JPMIN+1)
               FPTEMP= ZER
               DO 210 K=KF,KL
  210          FPTEMP= FPTEMP+FA(NBC+NLP(K),X12+ZLP(K))*CLP(K)
  220    FPSUM = FPSUM+FPTEMP*DUM
C
C NOW LOOP OVER THE ANGULAR FUNCTIONS, RETRIEVING THE ANGULAR
C INTEGRAL, MULTIPLYING BY THE RADIAL SUM AND STORE IN G
C
         FPSQ2 = FPSUM*Q2
         MMAX= JAMAX
         NN= 1
         DO 240 I=IAMIN,IAMAX
            DUMI= FPSQ2
            IF(NORMI) DUMI= DUMI*ECPFAC(I)
            IF(IANDJ) MMAX= I
            DO 240 J=JAMIN,MMAX
               DUMJ= DUMI
               IF(NORMJ) DUMJ= DUMJ*ECPFAC(J)
               SUM= ZER
               DO 230 M=-L,L
                  NLM=(NLM1-M)*NTLIM
  230          SUM= SUM+DCOEF4(NLM+I)*DCOEF4(NLM+J)
               G(NN)= G(NN)+DUMJ*SUM
  240    NN= NN+1
C
  300 CONTINUE
      RETURN
      END
C*MODULE ECP     *DECK ECPRA2
      SUBROUTINE ECPRA2(ICAB,NPNP,FP,COEFI,COEFJ,
     *                  DCOEF1,JFST1,LBECP1,DCOEF4,
     *                  DCOEF2,JFST2,LBECP2,G,
     *                  ZLM,LMF,LMX,LMY,LMZ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C ROUTINE ECPRA2 GATHERS ROUTINES ECPR12, ECPA12, ECPR22 AND ECPA22
C
      DIMENSION FP(*),COEFI(*),COEFJ(*),
     *          DCOEF1(*),JFST1(*),LBECP1(9,*),DCOEF4(*),
     *          DCOEF2(*),JFST2(*),LBECP2(6,*),G(*)
      DIMENSION ZLM(*),LMF(*),LMX(*),LMY(*),LMZ(*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
C ... ECP PARAMETERS ... GENERATED IN ECPPAR
C
      COMMON /ECP1  / X01,CAX,CAY,CAZ,CA,XCA,YCA,ZCA,
     *                X02,BAX,BAY,BAZ,BA,XBA,YBA,ZBA,
     *                PHASE,DAX,DAY,DAZ,DA,XDA,YDA,ZDA,XINT,KCNTR
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /ECP3  / ACX(7),ACY(7),ACZ(7),ABX(7),ABY(7),ABZ(7)
      COMMON /ECPDIM/ NCOEF1,NCOEF2,J1LEN,J2LEN,LLIM,NLIM,NTLIM,J4LEN
C IAMIN/MAX ARE THE ANGULAR   FUNCTION MIN AND MAX FOR I
C IPMIN/MAX ARE THE PRIMITIVE FUNCTION MIN AND MAX FOR I
C KF1/KL1 ARE THE FIRST AND LAST PRIMITIVES FORMING THE LMAX ECP
C     AND THUS THE TYPE 1 INTEGRAL
      LOGICAL IANDJ,NORM,NORMI,NORMJ
      COMMON /ECPIDX/ Q2,IAMIN,IAMAX,JAMIN,JAMAX,IPMIN,IPMAX,JPMIN,
     *                JPMAX,KF1,KL1,LLMX,NPC,NPB,IANDJ,NORM,NORMI,NORMJ
C
      COMMON /FICMN / ALF,XI,XP0,XP1
      COMMON /GBASE / NFST(8),NX(84),NY(84),NZ(84)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
C
      PARAMETER (ZER=0.0D+00, ONE=1.0D+00)
      PARAMETER (SR3=1.73205080756888D+00, SR5=2.23606797749979D+00)
      PARAMETER (SR7=2.64575131106459D+00)
      PARAMETER (S35=SR3*SR5, S57=SR5*SR7, S53=S57/SR3)
      PARAMETER (TOL=1.0D-10, FPI=12.566370614359D+00)
C
      DIMENSION ECPFAC(35),FIP(78),ZFNLM(121),CKL(11,11)
      DATA ECPFAC/7*ONE,3*SR3,3*ONE,6*SR5,1*S35,3*ONE,6*SR7,3*S53,3*S57/
C
C ROUTINE ECPR12
C
C TYPE 1 RADIAL INTEGRALS FOR THE GIVEN SET OF SHELLS
C SHELL SET TYPES <B|A|A> OR <A|A|C> (IE ICAB=(2 OR 3))
C     FP WILL CONTAIN THE RADIAL INTEGRALS UPON EXITING THIS ROUTINE
C     NPNP = N + N' INDICATES THE MAXIMUM VALUE FOR KAPPA NEEDED
C
C MAXIMUM NUMBER OF FP INTEGRALS NEEDED FOR THE GIVEN SET OF SHELLS
C FIP MUST BE DIMENSIONED LARGER THAN MAX NPNPMX (66 FOR G-HESS)
C
      NPNPMX=(NPNP*(NPNP+1))/2
C
C ZERO OUT THE INTEGRAL ARRAY
C
      CALL VCLR(FP,1,NPNPMX)
C
C LOOP OVER THE PRIMITIVES STORING THE RADIAL INTEGRALS IN FP
C
      DO 140 IG=IPMIN,IPMAX
         X01= EX(IG)
         DO 130 JG=JPMIN,JPMAX
            X02= EX(JG)
            X12= X01+X02
            DUM= COEFI(IG-IPMIN+1)*COEFJ(JG-JPMIN+1)
            IF(ICAB.EQ.2) THEN
               ALFA= X02*BA
               ALFX= ALFA*BA
            ELSE IF(ICAB.EQ.3) THEN
               ALFA= X01*CA
               ALFX= ALFA*CA
            END IF
            ALFI= ONE/(ALFA+ALFA)
            XP0 = EXP(-ALFX)
C
C LOOP OVER THE ECP PRIMITIVES
C
            DO 120 K=KF1,KL1
               XI = ONE/SQRT(X12+ZLP(K))
               ALF= ALFA*XI
               XP1= EXP(-ALFX+ALF*ALF)
               NLPK= NLP(K)
               CLPK= CLP(K)*DUM
C
C CALCULATE THE ACTUAL RADIAL INTEGRALS WITH FIECP
C
               CALL FIECP(FIP,ALFI,NLPK,NPNP-1)
C
C NOW STORE THEM IN FP AFTER MULTIPLYING BY THE COEFS OF EACH PRIMITIVE
C
               DO 110 N=1,NPNPMX
  110          FP(N)= FP(N)+FIP(N)*CLPK
  120       CONTINUE
C
  130    CONTINUE
  140 CONTINUE
C
C ROUTINE ECPA12
C
C ECP TYPE 1 ANGULAR INTEGRALS
C NOTE THE ACTUAL INTEGRAL HAS ALREADY BEEN DONE, WE JUST NEED
C TO COMPUTE THE APPROPRIATE FACTORS AND MULTIPLY BY THE STORED
C INTEGRAL.
C     FP WILL CONTAIN ALL THE REQUIRED RADIAL INTEGRALS UPON ENTRY
C     G WILL CONTAIN THE GRADIENT VECTOR UPON EXIT (NOT SYMMETRIZED)
C     LBECP1 CONTAINS THE PACKED INDEXES FOR THE ANGULAR INTEGRALS
C     DCOEF1 CONTAINS THE PACKED ANGULAR INTEGRAL VALUES
C     NPNP = N+N' GIVES THE MAXIMUM ANGULAR MOMENTUM (KA) NEEDED
C
      CALL ZFN(ZFNLM,NPNP-1,ZLM,LMF,LMX,LMY,LMZ)
C
C NN INDEXES THE GRADIENT VECTOR ARRAY
C
      FPIQ2 = FPI*Q2
      MMAX= JAMAX
      NN= 1
      DO 240 I=IAMIN,IAMAX
         DUMI= FPIQ2
         IF(NORMI) DUMI= DUMI*ECPFAC(I)
         IF(IANDJ) MMAX= I
         DO 240 J=JAMIN,MMAX
            DUMJ= DUMI
            IF(NORMJ) DUMJ= DUMJ*ECPFAC(J)
            IF(I.GE.J) THEN
               INDX= J+(I*(I-1))/2
            ELSE
               INDX= I+(J*(J-1))/2
            END IF
            JF = JFST1(INDX)
            JL = JFST1(INDX+1)-1
C
C ZERO OUT CKL WHICH STORES THE ANGULAR INTEGRALS TO BE USED LATER
C
            DO 210 L=1,NPNP
               DO 210 K=L,NPNP
  210       CKL(K,L)= ZER
C
C NOW CALCULATE THE ANGULAR INTEGRAL WHICH DOES NOT DEPEND ON PRIMITIVES
C
            NXYZ = NX(I)+NY(I)+NZ(I)+NX(J)+NY(J)+NZ(J)
            K= I
            M= 3
            IF(ICAB.EQ.2) THEN
               K= I+J-K
               M= 3+6-M
            END IF
            IF(I.GE.J) M= 3+6-M
            DO 220 JJ=JF,JL
               KA = LBECP1(1,JJ)+1
               LA = LBECP1(2,JJ)+1
               MU = LBECP1(3,JJ)
               KX = LBECP1(M+1,JJ)+NX(K)
               KY = LBECP1(M+2,JJ)+NY(K)
               KZ = LBECP1(M+3,JJ)+NZ(K)
               IF((KX+KY+KZ).NE.NXYZ) GO TO 220
               IX = LBECP1(4,JJ)+LBECP1(7,JJ)
               IY = LBECP1(5,JJ)+LBECP1(8,JJ)
               IZ = LBECP1(6,JJ)+LBECP1(9,JJ)
               LU = LA*(LA-1)-MU+1
               IF(ICAB.EQ.2) THEN
                  DT= ABX(KX+1-IX)*ABY(KY+1-IY)*ABZ(KZ+1-IZ)*ZFNLM(LU)
               ELSE
                  DT= ACX(KX+1-IX)*ACY(KY+1-IY)*ACZ(KZ+1-IZ)*ZFNLM(LU)
               END IF
               CKL(KA,LA)= CKL(KA,LA)+DCOEF1(JJ)*DT
  220       CONTINUE
C
C        -----  COMBINE THE RADIAL INTEGRALS WITH STRUCTURE     -----
C        -----  DEPENDENT COEFF.                                -----
C
C SUM IS USED TO COLLECT EACH TOTAL INTEGRAL
C
            SUM= ZER
            N= 1
            DO 230 K=1,NPNP
               DO 230 L=1,K
                  IF(ABS(CKL(K,L)).GT.TOL) SUM= SUM+FP(N)*CKL(K,L)
  230       N= N+1
C
C STORE THE COMPUTED INTEGRAL INTO G
C
            G(NN)= G(NN)+DUMJ*SUM
  240 NN= NN+1
C
      IF(LLMX.LT.2) RETURN
C
C ROUTINES ECPR22 AND ECPA22
C
      DO 500 LL=2,LLMX
C
C ROUTINE ECPR22
C
C TYPE 2 RADIAL INTEGRALS FOR THE GIVEN SET OF SHELLS
C SHELL SET TYPES <B|A|A> OR <A|A|C> (IE ICAB=(2 OR 3))
C     FP WILL CONTAIN THE RADIAL INTEGRALS UPON EXITING THIS ROUTINE
C     NPNP = N + N' INDICATES THE MAXIMUM VALUE FOR KA NEEDED
C     LL GIVES THE VALUE OF L FOR THE <|ULMAX-UL|> INTEGRAL
C
C ZERO OUT THE INTEGRAL ARRAY
C
         CALL VCLR(FP,1,NPNPMX)
C
C SET UP THE INDICES BASED ON LL
C
         KF = KFRST(KCNTR,LL)
         KL = KLAST(KCNTR,LL)
C
C LOOP OVER THE PRIMITIVES STORING THE RADIAL INTEGRALS IN FP
C
         DO 340 IG=IPMIN,IPMAX
            X01= EX(IG)
            DO 330 JG=JPMIN,JPMAX
               X02= EX(JG)
               X12= X01+X02
               DUM= COEFI(IG-IPMIN+1)*COEFJ(JG-JPMIN+1)
               IF(ICAB.EQ.2) THEN
                  GAMA= X02*BA
                  GAMX= GAMA*BA
               ELSE IF(ICAB.EQ.3) THEN
                  GAMA= X01*CA
                  GAMX= GAMA*CA
               END IF
               GAMI= ONE/(GAMA+GAMA)
               XP0 = EXP(-GAMX)
C LOOP OVER THE ECP PRIMITIVES
               DO 320 K=KF,KL
                  XI = ONE/SQRT(X12+ZLP(K))
                  ALF= GAMA*XI
                  XP1= EXP(-GAMX+ALF*ALF)
                  NLPK = NLP(K)
                  CLPK = CLP(K)*DUM
                  CALL FIECP(FIP,GAMI,NLPK,NPNP-1)
                  DO 310 N=1,NPNPMX
  310             FP(N)= FP(N)+FIP(N)*CLPK
  320          CONTINUE
  330       CONTINUE
  340    CONTINUE
C
C ROUTINE ECPA22
C
C ECP TYPE 2 ANGULAR INTEGRALS
C NOTE THE ACTUAL INTEGRAL HAS ALREADY BEEN DONE, WE JUST NEED
C TO COMPUTE THE APPROPRIATE FACTORS AND MULTIPLY BY THE STORED
C INTEGRAL.
C     FP WILL CONTAIN ALL THE REQUIRED RADIAL INTEGRALS UPON ENTRY
C     G WILL CONTAIN THE GRADIENT VECTOR UPON EXIT (NOT SYMMETRIZED)
C     LBECP2 CONTAINS THE PACKED INDEXES FOR THE ANGULAR INTEGRALS
C     DCOEF2 CONTAINS THE PACKED ANGULAR INTEGRAL VALUES
C     NPNP = N+N' GIVES THE MAXIMUM ANGULAR MOMENTUM (KA) NEEDED
C
         L2PL=(LL-2)*(LL-1)
         CALL ZFN(ZFNLM,NPNP-1,ZLM,LMF,LMX,LMY,LMZ)
C
C NN INDEXES THE GRADIENT VECTOR ARRAY
C
         FPIQ2 = FPI*Q2
         MMAX= JAMAX
         NN= 1
         DO 440 I=IAMIN,IAMAX
            DUMI= FPIQ2
            IF(NORMI) DUMI= DUMI*ECPFAC(I)
            IF(IANDJ) MMAX= I
            DO 440 J=JAMIN,MMAX
               DUMJ= DUMI
               IF(NORMJ) DUMJ= DUMJ*ECPFAC(J)
C
C ZERO OUT CKL WHICH STORES THE ANGULAR INTEGRALS TO BE USED LATER
C
               DO 400 L=1,NPNP
                  DO 400 K=L,NPNP
  400          CKL(K,L)= ZER
C
C NOW CALCULATE THE ANGULAR INTEGRAL WHICH DOES NOT DEPEND ON PRIMITIVES
C
               MX = NX(I)+NX(J)+1
               MY = NY(I)+NY(J)+1
               MZ = NZ(I)+NZ(J)+1
               IF(ICAB.EQ.2) THEN
                  K= I
               ELSE
                  K= J
               END IF
               KT = NX(K)+NY(K)+NZ(K)+1
               L= LL-2
               DO 420 M=L,-L,-1
                  NLM=(L2PL-M)*NTLIM
                  CT = DCOEF4(NLM+K)
                  IF(ABS(CT).LT.TOL) GO TO 420
                  NLM= NLM+I+J
                  JF = JFST2(NLM-K)
                  JL = JFST2(NLM-K+1)-1
                  DO 410 JJ=JF,JL
                     LA = LBECP2(1,JJ)+1
                     KA = LBECP2(2,JJ)+KT
                     MU = LBECP2(3,JJ)
                     KX = LBECP2(4,JJ)+NX(K)
                     KY = LBECP2(5,JJ)+NY(K)
                     KZ = LBECP2(6,JJ)+NZ(K)
                     LU = LA*(LA-1)-MU+1
                     IF(ICAB.EQ.2) THEN
                        DT=CT*ABX(MX-KX)*ABY(MY-KY)*ABZ(MZ-KZ)*ZFNLM(LU)
                     ELSE
                        DT=CT*ACX(MX-KX)*ACY(MY-KY)*ACZ(MZ-KZ)*ZFNLM(LU)
                     END IF
                     CKL(KA,LA)= CKL(KA,LA)+DCOEF2(JJ)*DT
  410             CONTINUE
  420          CONTINUE
C
C COMBINE THE ANGULAR AND THE RADIAL INTEGRALS
C
               SUM= ZER
               N= 1
               DO 430 K=1,NPNP
                  DO 430 L=1,K
               IF(ABS(CKL(K,L)).GT.TOL) SUM= SUM+FP(N)*CKL(K,L)
  430          N= N+1
C
C STORE THE COMPUTED INTEGRAL INTO G
C
               G(NN)= G(NN)+DUMJ*SUM
  440    NN= NN+1
C
  500 CONTINUE
C
      RETURN
      END
C*MODULE ECP     *DECK ECPDRA
      SUBROUTINE ECPDRA(IC4C,NPNP,FP,FQ,COEFI,COEFQ,COEFJ,
     *                  DCOEF1,JFST1,LBECP1,
     *                  DCOEF2,JFST2,LBECP2,G,
     *                  ZLM,LMF,LMX,LMY,LMZ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C ECPDRA GATHERS ROUTINES ECPD14, ECPR14, ECPA14, ECPR24 AND ECPA24
C
      DIMENSION FP(*),FQ(*),COEFI(*),COEFQ(*),COEFJ(*),
     *          DCOEF1(*),JFST1(*),LBECP1(9,*),
     *          DCOEF2(*),JFST2(*),LBECP2(6,*),G(*)
      DIMENSION ZLM(*),LMF(*),LMX(*),LMY(*),LMZ(*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
C ... ECP PARAMETERS ... GENERATED IN ECPPAR
C
      COMMON /ECP1  / X01,CAX,CAY,CAZ,CA,XCA,YCA,ZCA,
     *                X02,BAX,BAY,BAZ,BA,XBA,YBA,ZBA,
     *                PHASE,DAX,DAY,DAZ,DA,XDA,YDA,ZDA,XINT,KCNTR
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /ECP3  / ACX(7),ACY(7),ACZ(7),ABX(7),ABY(7),ABZ(7)
      LOGICAL CANDB
      COMMON /ECP4  / P12(3,2),R12,ACO(3),CANDB
      COMMON /ECPDIM/ NCOEF1,NCOEF2,J1LEN,J2LEN,LLIM,NLIM,NTLIM,J4LEN
C IAMIN/MAX ARE THE ANGULAR   FUNCTION MIN AND MAX FOR I
C IPMIN/MAX ARE THE PRIMITIVE FUNCTION MIN AND MAX FOR I
C KF1/KL1 ARE THE FIRST AND LAST PRIMITIVES FORMING THE LMAX ECP
C     AND THUS THE TYPE 1 INTEGRAL
      LOGICAL IANDJ,NORM,NORMI,NORMJ
      COMMON /ECPIDX/ Q2,IAMIN,IAMAX,JAMIN,JAMAX,IPMIN,IPMAX,JPMIN,
     *                JPMAX,KF1,KL1,LLMX,NPC,NPB,IANDJ,NORM,NORMI,NORMJ
C ECP COMMON STUFF
      COMMON /ZFNCM / X,Y,Z
C
      COMMON /FICMN / ALF,XI,XP0,XP1
      COMMON /FJCMN / ALEF,BEIT,XXI,XPLS,XMNS,XP
      COMMON /FJNEW / XKA,XKB,GAM1,GAM2,A1,A2,C
      COMMON /GBASE / NFST(8),NX(84),NY(84),NZ(84)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
C
      DIMENSION FIP(78),ZFNLM(121),CKLU(23,12,12)
      DIMENSION CKLC(11,11),CKLB(11,11),CPQ(11,11,11)
      DIMENSION FJPQ(11,11,11),ZFNLMC(121),ZFNLMB(121)
      SAVE FJPQ,ZFNLMC,ZFNLMB
C
      PARAMETER (ZER=0.0D+00, ONE=1.0D+00, ABLIM=1.0D-01)
      PARAMETER (SR3=1.73205080756888D+00, SR5=2.23606797749979D+00)
      PARAMETER (SR7=2.64575131106459D+00)
      PARAMETER (S35=SR3*SR5, S57=SR5*SR7, S53=S57/SR3)
      PARAMETER (TOL=1.0D-10, FPI=12.566370614359D+00)
      PARAMETER (TM6=1.0D-06, ONDS4P=0.28209479177388D+00)
C
      DIMENSION ECPFAC(35)
      DATA ECPFAC/7*ONE,3*SR3,3*ONE,6*SR5,1*S35,3*ONE,6*SR7,3*S53,3*S57/
C
      IF(IC4C.GT.0) GO TO 200
C
C ROUTINE ECPD14, FOR IC4C < 0
C ROUTINE ECPR14, FOR IC4C = 0
C
C TYPE 1 RADIAL INTEGRALS FOR THE GIVEN SET OF SHELLS
C SHELL SET TYPES <B|A|C> (IE ICAB=4)
C
C MAXIMUM NUMBER OF FP INTEGRALS NEEDED FOR THE GIVEN SET OF SHELLS
C
      NPNPMX=(NPNP*(NPNP+1)*(2*NPNP+1))/6
C
C ZERO OUT THE INTEGRAL ARRAY
C
      CALL VCLR(FP,1,NPNPMX)
      IF(IC4C.LT.0)
     *CALL VCLR(FQ,1,NPNPMX)
C
C LOOP OVER THE PRIMITIVES STORING THE RADIAL INTEGRALS IN FP
C
      DO 160 IG=IPMIN,IPMAX
         X01= EX(IG)
         DO 160 JG=JPMIN,JPMAX
            X02= EX(JG)
            X12= X01+X02
            DUM= COEFI(IG-IPMIN+1)*COEFJ(JG-JPMIN+1)
            IF(IC4C.LT.0)
     *      DUQ= COEFQ(IG-IPMIN+1)*COEFJ(JG-JPMIN+1)
C
C TRANSFORM THE GAUSSIANS TO CENTER D ACTUALLY C AND B MAY BE THE
C SAME AS LONG AS NEITHER IS THE SAME AS A. IF C AND B ARE DIFFERENT
C THEN A PHASE FACTOR IS ALSO REQUIRED
C
            X21= ONE/X12
            Y01= X01*X21
            Y02= ONE-Y01
            Y12= Y01*X02
            DAX= P12(1,1)+(P12(1,2)-P12(1,1))*Y02-ACO(1)
            DAY= P12(2,1)+(P12(2,2)-P12(2,1))*Y02-ACO(2)
            DAZ= P12(3,1)+(P12(3,2)-P12(3,1))*Y02-ACO(3)
            DA = SQRT(DAX*DAX+DAY*DAY+DAZ*DAZ)
            IF(.NOT.CANDB) THEN
               DUM= DUM*EXP(-R12*Y12)
               IF(IC4C.LT.0)
     *         DUQ= DUQ*EXP(-R12*Y12)
            END IF
C
C IF CENTER D AND A BOTH LIE ON THE ORIGIN THEN ALL INTEGRALS FOR
C LA<>0 ARE 0 SO SPECIAL CASE THAT CODE
C
            IF(DA.GE.TM6) THEN
               ALFA= X12*DA
               ALFX= ALFA*DA
               ALFI= ONE/(ALFA+ALFA)
               X  = DAX/DA
               Y  = DAY/DA
               Z  = DAZ/DA
               CALL ZFN(ZFNLM,NPNP-1,ZLM,LMF,LMX,LMY,LMZ)
               XP0 = EXP(-ALFX)
               DO 130 K=KF1,KL1
                  XI  = ONE/SQRT(X12+ZLP(K))
                  ALF = ALFA*XI
                  XP1 = EXP(-ALFX+ALF*ALF)
                  NLPK= NLP(K)
                  CLPK= CLP(K)*DUM
                  IF(IC4C.LT.0)
     *            CLPQ= CLP(K)*DUQ
                  CALL FIECP(FIP,ALFI,NLPK,NPNP-1)
                  N= 1
                  NN= 1
                  DO 120 KK=1,NPNP
                     DO 120 L=1,KK
                        FIPTEM= FIP(NN)*CLPK
                        IF(IC4C.LT.0)
     *                  FIPTEQ= FIP(NN)*CLPQ
                        KKLL= L*(L-1)+1+L
                        DO 110 MU=1,2*L-1
                           FP(N)= FP(N)+FIPTEM*ZFNLM(KKLL-MU)
                           IF(IC4C.LT.0)
     *                     FQ(N)= FQ(N)+FIPTEQ*ZFNLM(KKLL-MU)
  110                   N= N+1
  120             NN= NN+1
  130          CONTINUE
            ELSE
C SPECIAL CASE WHEN D=A=0. THERE ARE ONLY L=0 INTEGRALS, BUT
C THEY MUST BE STORED THE SAME AS ABOVE SO INDEX THEM INTO FP.
C THEY ALSO MUST BE MULTIPLIED BY THE 0 SPHERICAL HAMONIC 1/SQRT(4PI)
               DO 150 K=KF1,KL1
                  ZETA= X12+ZLP(K)
                  NLPK= NLP(K)
                  CLPK= CLP(K)*ONDS4P*DUM
                  IF(IC4C.LT.0)
     *            CLPQ= CLP(K)*ONDS4P*DUQ
                  DO 140 N=1,NPNP
                     NK=(N*(N-1)*(2*N-1))/6+1
                     FATEMP= FA(N+NLPK-1,ZETA)
                     FP(NK)= FP(NK)+FATEMP*CLPK
                     IF(IC4C.LT.0)
     *               FQ(NK)= FQ(NK)+FATEMP*CLPQ
  140             CONTINUE
  150          CONTINUE
            END IF
  160 CONTINUE
C
  200 CONTINUE
C
C ROUTINE ECPA14
C
C ECP TYPE 1 ANGULAR INTEGRALS
C ICAB = 4 GENERAL CASE WITH 3 DIFFERENT CENTERS <CAB>
C NOTE THE ACTUAL INTEGRAL HAS ALREADY BEEN DONE, WE JUST NEED
C TO COMPUTE THE APPROPRIATE FACTORS AND MULTIPLY BY THE STORED
C INTEGRAL.
C     FP WILL CONTAIN ALL THE REQUIRED RADIAL INTEGRALS UPON ENTRY
C     G WILL CONTAIN THE GRADIENT VECTOR UPON EXIT (NOT SYMMETRIZED)
C     LBECP1 CONTAINS THE PACKED INDEXES FOR THE ANGULAR INTEGRALS
C     DCOEF1 CONTAINS THE PACKED ANGULAR INTEGRAL VALUES
C     NPNP = N+N' GIVES THE MAXIMUM ANGULAR MOMENTUM (KA) NEEDED
C
C NN INDEXES THE GRADIENT VECTOR ARRAY
C
      FPIQ2 = FPI*Q2
      MMAX= JAMAX
      NN= 1
      DO 245 I=IAMIN,IAMAX
         DUMI= FPIQ2
         IF(NORMI) DUMI= DUMI*ECPFAC(I)
         IF(IANDJ) MMAX= I
         DO 240 J=JAMIN,MMAX
            DUMJ= DUMI
            IF(NORMJ) DUMJ= DUMJ*ECPFAC(J)
            IF(I.GE.J) THEN
               INDX= J+(I*(I-1))/2
            ELSE
               INDX= I+(J*(J-1))/2
            END IF
            JF = JFST1(INDX)
            JL = JFST1(INDX+1)-1
            NXC= NX(I)+1
            NYC= NY(I)+1
            NZC= NZ(I)+1
            NXB= NX(J)+1
            NYB= NY(J)+1
            NZB= NZ(J)+1
C
C ZERO OUT CKL WHICH STORES THE ANGULAR INTEGRALS TO BE USED LATER
C THESE ZEROING LOOPS ORIGINALLY MATCHED THE COMPUTE LOOPS 150/140/130
C BELOW.  HOWEVER, THE VALUE OF LA+MU IN LOOP 120 CAN BE AS BIG
C AS NPNP+1, SO TO AVOID CORE DUMPS, WE ZERO MORE OF THE ARRAY OUT.
C ORIGINAL DIMENSION OF CKLU WAS (22,11,11), CHANGED TO (23,12,12).
C
            DO 210 K=1,NPNP+1
               DO 210 L=1,K
                  DO 210 MU=1,2*L-1
  210       CKLU(MU,L,K)= ZER
C
            DO 220 K=JF,JL
               KA = LBECP1(1,K)+1
               LA = LBECP1(2,K)+1
               MU = LBECP1(3,K)
               IF(I.GE.J) THEN
                  KX = LBECP1(4,K)
                  KY = LBECP1(5,K)
                  KZ = LBECP1(6,K)
                  KXP= LBECP1(7,K)
                  KYP= LBECP1(8,K)
                  KZP= LBECP1(9,K)
               ELSE
                  KXP= LBECP1(4,K)
                  KYP= LBECP1(5,K)
                  KZP= LBECP1(6,K)
                  KX = LBECP1(7,K)
                  KY = LBECP1(8,K)
                  KZ = LBECP1(9,K)
               END IF
      CKLU(LA+MU,LA,KA)= CKLU(LA+MU,LA,KA)+DCOEF1(K)*
     *                   ACX(NXC-KX) *ACY(NYC-KY) *ACZ(NZC-KZ) *
     *                   ABX(NXB-KXP)*ABY(NYB-KYP)*ABZ(NZB-KZP)
  220       CONTINUE
C
C SUM IS USED TO COLLECT EACH TOTAL INTEGRAL
C
            SUM= ZER
            N= 1
            DO 232 K=1,NPNP
               DO 231 L=1,K
                  DO 230 MU=1,2*L-1
                     CKLTEM= CKLU(MU,L,K)
                     IF(IC4C.LE.0) THEN
                        IF(ABS(CKLTEM).GT.TOL) SUM= SUM+FP(N)*CKLTEM
                     ELSE
                        IF(ABS(CKLTEM).GT.TOL) SUM= SUM+FQ(N)*CKLTEM
                     END IF
                     N= N+1
  230             CONTINUE
  231          CONTINUE
  232       CONTINUE
C
C STORE THE COMPUTED INTEGRAL INTO G
C
            G(NN)= G(NN)+DUMJ*SUM
            NN= NN+1
  240    CONTINUE
  245 CONTINUE
C
      IF(LLMX.LT.2) RETURN
C
C ROUTINES ECPR24 AND ECPA24
C
      FPISQ2= FPI*FPIQ2
      DO 500 LL=2,LLMX
         L= LL-2
C
C ROUTINE ECPR24
C
C TYPE 2 RADIAL INTEGRALS FOR THE GIVEN SET OF SHELLS
C SHELL SET TYPES <B|A|C> (IE ICAB=4)
C
C SETUP VARIOUS PARAMETERS FOR LOOPS
C
         NP1  = MAX0(NPC,NPB)
         NP1PL= NP1+L
         NPCPL= NPC+L
         NPBPL= NPB+L
         NMAX =(NPNP-1)*NPCPL*(NPBPL/2+1)+6
         LTMAX= MAX0(NPNP,NP1PL)
         LEMX = MAX0(L,LTMAX/2)
C
C SET UP THE INDICES BASED ON LL
C
         KF = KFRST(KCNTR,LL)
         KL = KLAST(KCNTR,LL)
C
C ZERO OUT THE INTEGRAL ARRAY
C
         CALL VCLR(FP,1,NMAX)
C
C LOOP OVER THE PRIMITIVES STORING THE RADIAL INTEGRALS IN FP
C
         DO 340 IG=IPMIN,IPMAX
            X01= EX(IG)
            ALFA= X01*CA
            ALFX= ALFA*CA
            ALFI= ONE/(ALFA+ALFA)
            DO 340 JG=JPMIN,JPMAX
               X02= EX(JG)
               X12= X01+X02
               BETA= X02*BA
               BETX= BETA*BA
               BETI= ONE/(BETA+BETA)
               ALBE= ALFX+BETX
               XP = EXP(-ALBE)
               DUM= COEFI(IG-IPMIN+1)*COEFJ(JG-JPMIN+1)
               DO 330 K=KF,KL
                  XXI= ONE/SQRT(X12+ZLP(K))
                  ALEF= ALFA*XXI
                  BEIT= BETA*XXI
                  IF(ALEF*BEIT.GT.ABLIM) THEN
                     XKA= ALFA+ALFA
                     XKB= BETA+BETA
                     A1 = ALEF+BEIT
                     A2 = ALEF-BEIT
                     C  = X12+ZLP(K)
                     XPLS= EXP(-ALBE+(ALEF+BEIT)**2)
                     XMNS= EXP(-ALBE+(ALEF-BEIT)**2)
                     GAM1= XPLS*0.25D+00
                     GAM2= XMNS*0.25D+00
                  ELSE
                     XPLS= EXP(-ALBE+ALEF*ALEF)
                     XMNS= EXP(-ALBE+BEIT*BEIT)
                  END IF
                  NLPK = NLP(K)
                  CLPK = CLP(K)*DUM
                  CALL FJECP(FJPQ,ALFI,BETI,NLPK,NPNP,LTMAX,LEMX)
C
C THE FIRST BLOCK IS SPECIAL SINCE THERE ARE ONLY VALUES ON THE
C DIAGONAL SO STICK INTO LINEAR MEMORY (MAX OF 6 ELEMENTS)
C
                  DO 310 IN=1,LEMX+1
                     FP(IN)= FP(IN)+FJPQ(IN,IN,1)*CLPK
  310             CONTINUE
C
C DO HIGHER BLOCKS HERE, NOTE: EVERY OTHER ELEMENT IS ZERO AND IS
C THEREFORE ZERO (ACTUALLY THE ANGULAR TERM WILL BE ZERO)
C
                  IF(NPNP.LT.2) GO TO 330
                  NBEG= 1
                  N= 6+1
                  DO 320 IN=2,NPNP
                     NBEG= 1-NBEG
                     LBEG= NBEG
                     DO 320 IP=1,NPCPL
                        LBEG= 1-LBEG
                        DO 320 IQ=1+LBEG,NPBPL,2
                  FP(N)= FP(N)+FJPQ(IQ,IP,IN)*CLPK
  320             N= N+1
  330          CONTINUE
  340    CONTINUE
C
C ROUTINE ECPA24
C
C ECP TYPE 2 ANGULAR INTEGRALS
C NOTE THE ACTUAL INTEGRAL HAS ALREADY BEEN DONE, WE JUST NEED
C TO COMPUTE THE APPROPRIATE FACTORS AND MULTIPLY BY THE STORED
C INTEGRAL.
C     FP WILL CONTAIN ALL THE REQUIRED RADIAL INTEGRALS UPON ENTRY
C     G WILL CONTAIN THE GRADIENT VECTOR UPON EXIT (NOT SYMMETRIZED)
C     LBECP2 CONTAINS THE PACKED INDEXES FOR THE ANGULAR INTEGRALS
C     DCOEF2 CONTAINS THE PACKED ANGULAR INTEGRAL VALUES
C     NPNP = N+N' GIVES THE MAXIMUM ANGULAR MOMENTUM (KA) NEEDED
C     L = THE ANGULAR MOMENTUM OF THE CORE REPRESENTED BY THE ECP
C
         L= LL-2
         L1MAX= MAX0(1,L+1)
         L2PL=(LL-1)*(LL-2)
         X = XCA
         Y = YCA
         Z = ZCA
         CALL ZFN(ZFNLMC,NPCPL-1,ZLM,LMF,LMX,LMY,LMZ)
         X = XBA
         Y = YBA
         Z = ZBA
         CALL ZFN(ZFNLMB,NPBPL-1,ZLM,LMF,LMX,LMY,LMZ)
C
C NN INDEXES THE INTEGRAL VECTOR ARRAY
C
CC       FPISQ2= FPI*FPI*Q2
         MMAX= JAMAX
         NN= 1
         DO 445 I=IAMIN,IAMAX
            DUMI = FPISQ2
            IF(NORMI) DUMI= DUMI*ECPFAC(I)
            IF(IANDJ) MMAX= I
            DO 440 J=JAMIN,MMAX
               DUMJ= DUMI
               IF(NORMJ) DUMJ= DUMJ*ECPFAC(J)
               NXC = NX(I)+1
               NYC = NY(I)+1
               NZC = NZ(I)+1
               NXB = NX(J)+1
               NYB = NY(J)+1
               NZB = NZ(J)+1
C
C ZERO OUT CKL WHICH STORES THE ANGULAR INTEGRALS TO BE USED LATER
C
               DO 410 M=1,2*NP1-1
                  DO 410 N=1,NPCPL
                     DO 410 K=1,NPBPL
  410          CPQ(K,N,M)= ZER
C
C NOW CALCULATE THE ANGULAR INTEGRAL WHICH DOES NOT DEPEND ON PRIMITIVES
C
               NLM=(L2PL-L)*NTLIM
               DO 420 M=L,-L,-1
                  DO 412 N=1,NP1
                     DO 412 K=1,NP1PL
                        CKLC(K,N)= ZER
                        CKLB(K,N)= ZER
  412             CONTINUE
                  JF = JFST2(NLM+I)
                  JL = JFST2(NLM+I+1)-1
CBB               IF(NLM+I.GT.J2LEN .OR. JL.GT.NCOEF2) THEN
CBB                  WRITE(6,*) 'NCOEF2 OUT A24',NLM,I,NCOEF2,JF,JL
CBB                  CALL ABRT
CBB               END IF
                  DO 414 JJ=JF,JL
                     LA = LBECP2(1,JJ)+1
                     KA = LBECP2(2,JJ)+1
                     MU = LBECP2(3,JJ)
                     KX = LBECP2(4,JJ)
                     KY = LBECP2(5,JJ)
                     KZ = LBECP2(6,JJ)
                     LU = LA*(LA-1)-MU+1
                     CKLC(LA,KA)= CKLC(LA,KA)+DCOEF2(JJ)*ZFNLMC(LU)*
     *                            ACX(NXC-KX)*ACY(NYC-KY)*ACZ(NZC-KZ)
  414             CONTINUE
                  JF = JFST2(NLM+J)
                  JL = JFST2(NLM+J+1)-1
CBB               IF(NLM+J.GT.J2LEN .OR. JL.GT.NCOEF2) THEN
CBB                  WRITE(6,*) 'NCOEF2 OUT A24',NLM,J,NCOEF2,JF,JL
CBB                  CALL ABRT
CBB               END IF
                  DO 416 JJ=JF,JL
                     LA = LBECP2(1,JJ)+1
                     KA = LBECP2(2,JJ)+1
                     MU = LBECP2(3,JJ)
                     KX = LBECP2(4,JJ)
                     KY = LBECP2(5,JJ)
                     KZ = LBECP2(6,JJ)
                     LU = LA*(LA-1)-MU+1
                     CKLB(LA,KA)= CKLB(LA,KA)+DCOEF2(JJ)*ZFNLMB(LU)*
     *                            ABX(NXB-KX)*ABY(NYB-KY)*ABZ(NZB-KZ)
  416             CONTINUE
C        -----  MULTIPLY COEFF OF THE TWO PROJECTION HALVES     -----
C        -----  AND STORE.                                      -----
                  DO 418 K=1,NP1
                     DO 418 N=1,NPCPL
                        CT= CKLC(N,K)
                        DO 418 K2=K,K+NP1-1
                           DO 418 N2=1,NPBPL
  418             CPQ(N2,N,K2)= CPQ(N2,N,K2)+CT*CKLB(N2,K2+1-K)
  420          NLM= NLM+NTLIM
C
C        -----  COMBINE THE RADIAL INTEGRALS WITH STRUCTURE     -----
C        -----  DEPENDENT COEFF.                                -----
C
C DO KA=0 FIRST SINCE IT ONLY HAS VALUES ON THE DIAGONAL
C SUM IS USED TO COLLECT EACH TOTAL INTEGRAL
C
               SUM= ZER
               DO 422 N=1,L1MAX
                  SUM= SUM+FP(N)*CPQ(N,N,1)
  422          CONTINUE
C
C IF THERE ARE KA VALUES LARGER THAN 0 DO THOSE NOW
C NOTE: EVERY OTHER ANGULAR INTEGRAL WILL BE ZERO BY SYMMETRY AND IS
C THEREFORE SKIPPED
C
               IF(NPNP.LT.2) GO TO 425
               JBEG= 1
               N= 6+1
               DO 424 KK=2,NPNP
                  JBEG= 1-JBEG
                  LBEG= JBEG
                  DO 424 JJ=1,NPCPL
                     LBEG= 1-LBEG
                     DO 424 M=LBEG+1,NPBPL,2
               IF(ABS(CPQ(M,JJ,KK)).GT.TOL) SUM= SUM+FP(N)*CPQ(M,JJ,KK)
  424          N= N+1
  425          CONTINUE
C
C STORE THE COMPUTED INTEGRAL INTO G
C
               G(NN)= G(NN)+DUMJ*SUM
               NN= NN+1
  440       CONTINUE
  445    CONTINUE
C
  500 CONTINUE
C
      RETURN
      END
C*MODULE ECP     *DECK ECPINT
      SUBROUTINE ECPINT(HECP,H0,DCOEF1,JFST1,LBECP1,DCOEF4,DCOEF2,
     *                  JFST2,LBECP2,FPQR,ZLM,LMF,LMX,LMY,LMZ,L2,INAT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION HECP(L2),H0(L2),DCOEF1(*),JFST1(*),LBECP1(9,*),
     *          DCOEF4(*),DCOEF2(*),JFST2(*),LBECP2(6,*),
     *          FPQR(25,25,25),ZLM(*),LMF(*),LMX(*),LMY(*),LMZ(*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXATM=2000,
     *           MXAO=8192)
C
C ... ECP PARAMETERS ... GENERATED IN ECPPAR
C
      COMMON /ECP1  / X01,CAX,CAY,CAZ,CA,XCA,YCA,ZCA,
     *                X02,BAX,BAY,BAZ,BA,XBA,YBA,ZBA,
     *                PHASE,DAX,DAY,DAZ,DA,XDA,YDA,ZDA,XINT,KCNTR
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      LOGICAL CANDB
      COMMON /ECP4  / P12(3,2),R12,ACO(3),CANDB
      COMMON /ECPDIM/ NCOEF1,NCOEF2,J1LEN,J2LEN,LLIM,NLIM,NTLIM,J4LEN
C IAMIN/MAX ARE THE ANGULAR   FUNCTION MIN AND MAX FOR I
C IPMIN/MAX ARE THE PRIMITIVE FUNCTION MIN AND MAX FOR I
C KF1/KL1 ARE THE FIRST AND LAST PRIMITIVES FORMING THE LMAX ECP
C     AND THUS THE TYPE 1 INTEGRAL
      LOGICAL IANDJ,NORM,NORMI,NORMJ
      COMMON /ECPIDX/ Q2,IAMIN,IAMAX,JAMIN,JAMAX,IPMIN,IPMAX,JPMIN,
     *                JPMAX,KF1,KL1,LLMX,NPC,NPB,IANDJ,NORM,NORMI,NORMJ
C
      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 /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      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
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
C THE FOLLOWING ARE TO MAKE USE OF SYMMETRY
C
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
C FP IS USED FOR RADIAL INTEGRAL STORAGE, MUST BE AT LEAST 2*11**3
C G-INTEGRAL STORAGE - MUST BE IMAX*JMAX(15*28 FOR G-HESS)
C 7 - FOR DIMENSION OF NLIM (FOR G-HESS)
C 13 - DIMENSIONED UP TO NMAX=2*NLIM-1
C NEED NX ETC DIMENSIONED UP TO 35 FOR G ENERGY 74 FOR G HESS
C
      DIMENSION FP(2*11*11*11),G(28*15),COEFI(MXGSH),COEFJ(MXGSH)
C
C DUMMY ARGUMENTS FOR ECPDRA WHEN IC4C=0
      DIMENSION FQ(1),COEFQ(1)
C
      PARAMETER (ZER=0.0D+00)
C
      LOGICAL CANDA,AANDB,DBUG,NXT
      DOUBLE PRECISION INT1,MOROKM
C
C IANG GIVES THE MAXIMUM ANGULAR MOMENTUM FOR EACH SHELL
C
      DIMENSION IANG(35),MI(48)
      DATA IANG/1,3*2,6*3,10*4,15*5/
      DATA ECPIN/8HECPINT  /, INT1/8HINT1     /, MOROKM/8HMOROKUMA/
C
C     ----- ECP MODIFICATIONS TO H INTEGRALS -----
C
C        -----  ROUTINE CONTROLS THE CALCULATION OF THE LOCAL    -----
C        -----  POTENTIAL.  -ECP- PROGRAM  BY L. KAHN.           -----
C        -----  GAUSSIAN BASIS INTEGRALS.                        -----
C        -----  WRITTEN AT BATTELLE, OCT. 1972.                  -----
C        -----  UPDATE IN  1978.                                 -----
C        -----  REWRITTEN AT IMS, OKAZAKI JAPAN  1980.           -----
C        -----  GRADSCF VERSION AUGUST 1981.                     -----
C        -----  HONDO VERSION APRIL 1984.                        -----
C        -----  GAMESS VERSION OCTOBER 1989.                     -----
C        -----  REWRITTEN FOR GAMESS JUNE 1995 BMB               -----
C
      FQ(1)   =ZER
      COEFQ(1)=ZER
C
C     ----- INITIALIZE PARALLEL -----
C
      NXT = IBTYP.EQ.1
      IPCOUNT= ME-1
      NEXT=-1
      LCNT=-1
      IF(GOPARR) CALL VCLR(HECP,1,L2)
C
      DBUG= EXETYP.EQ.ECPIN .OR. EXETYP.EQ.INT1
      DBUG= MASWRK .AND.(NPRINT.EQ.3 .OR. DBUG)
C
C ALWAYS WANT SHELLS NORMALIZED WITH REGULAR INTEGRALS
C
      NORM = NORMF.NE.1 .OR. NORMP.NE.1
      NORMI= NORM
      NORMJ= NORM
      J= 0
      DO 110 I=1,NUM
C        IA(I)=(I*(I-1))/2
         IA(I)= J
  110 J= J+I
C
C     -----  READ IN ECP FORMULAS AND DATA  -----
C
      LDAF91=(J1LEN-1)/NWDVAR+1+(9*NCOEF1-1)/NWDVAR+1
      LDAF93=(J2LEN-1)/NWDVAR+1+(6*NCOEF2)/NWDVAR
      CALL DAREAD(IDAF,IODA,DCOEF1,NCOEF1,90,0)
      CALL DAREAD(IDAF,IODA,JFST1 ,LDAF91,91,1)
      CALL DAREAD(IDAF,IODA,DCOEF2,NCOEF2,92,0)
      CALL DAREAD(IDAF,IODA,JFST2 ,LDAF93,93,1)
      CALL DAWT
      CALL ERRT
      CALL DAWERT
      CALL ECPINI(LMF,LMX,LMY,LMZ)
      CALL ZTAB(ZLM)
      CALL FTAB(FPQR,NLIM-1)
      CALL ECCOD3(FPQR,DCOEF4,ZLM,LMF,LMX,LMY,LMZ)
      IBEG = 1
      IEND = NSHELL
      JBEG = 1
      LOCIJ= 0
      NATST= 1
      NATED= NAT
      ISAVE= 0
C
C INAT WILL BE NON-ZERO FOR LCD PARTS OF A RUN
C WHICH REQUIRE NO SYMMETRY
C
      IF(INAT.EQ.0) THEN
C
C     ----- SET PARAMETERS FOR NORMAL OR MOROKUMA TYPE RUNS -----
C           ISAVE .EQ. 0 : SAVE H+HECP TO DAF 11
C           ISAVE .EQ. 1 : SAVE H+HECP TO DAF 11 AND TO DAF 311
C
         IF(RUNTYP.EQ.MOROKM)
     *   CALL STINT1(IBEG,IEND,JBEG,LOCIJ,NATST,NATED,ISAVE,L1X,L2X)
         LL2 =(NUM*NUM+NUM)/2
      ELSE
         NATST= INAT
         NATED= INAT
      END IF
C
C INITIALIZE HECP---NEEDED WHEN SYMMETRY IS USED
C
      NDUM=(NUM*(NUM+1))/2
      CALL VCLR(HECP,1,NDUM)
C
C        -----  ISHELL  -----
C
      DO 240 II=IBEG,IEND
         IF(INAT.EQ.0) THEN
C CHECK TO SEE IF SYMMETRY CAN BE USED TO AVOID THIS SHELL
            DO 120 IT=1,NT
               ID=MAPSHL(II,IT)
               IF(ID.GT.II) GO TO 240
               MI(IT)=ID
  120       CONTINUE
         END IF
C
C        -----  GO PARALLEL! -----
C
         IF(NXT .AND. GOPARR) THEN
            LCNT= LCNT+1
            IF(LCNT.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
            IF(LCNT.NE.NEXT) GO TO 240
         END IF
C
         I1 = KSTART(II)
         I2 = I1+KNG(II)-1
         IPMIN = I1
         IPMAX = I2
         ICNTR = KATOM(II)
         IMIN = KMIN(II)
         IMAX = KMAX(II)
         LOCI = KLOC(II)-IMIN-LOCIJ
         IAMIN = IMIN
         IAMAX = IMAX
         IIMAX = 1
         IF(IMIN.EQ.1 .AND. IMAX.EQ.4) IIMAX = 2
         DO 230 III=1,IIMAX
            IF(IIMAX.EQ.2) THEN
               IF(III.EQ.1) THEN
                  IAMIN = 1
                  IAMAX = 1
               ELSE
                  IAMIN = 2
                  IAMAX = 4
               END IF
            END IF
            NPC0= IANG(IAMAX)
            DO 130 IG=IPMIN,IPMAX
               IF(IAMIN.LE.35) T01= CG(IG)
               IF(IAMIN.LE.20) T01= CF(IG)
               IF(IAMIN.LE.10) T01= CD(IG)
               IF(IAMIN.LE. 4) T01= CP(IG)
               IF(IAMIN.EQ. 1) T01= CS(IG)
               COEFI(IG-IPMIN+1)= T01
  130       CONTINUE
C
C        -----  JSHELL  -----
C
            DO 200 JJ=JBEG,II
               Q2=1
               IF(INAT.EQ.0) THEN
C CHECK SYMMETRY AGAIN
                  N2=0
                  DO 140 IT=1,NT
                     JD=MAPSHL(JJ,IT)
                     IF(JD.GT.II) GO TO 200
                     ID=MAX0(MI(IT),JD)
                     JD=MIN0(MI(IT),JD)
                     IF(ID.LT.II) GO TO 140
                     IF(JD.LT.JJ) GO TO 140
                     IF(JD.GT.JJ) GO TO 200
                     N2=N2+1
  140             CONTINUE
                  Q2 = NT
                  Q2 = Q2/N2
               END IF
C
C        -----  GO PARALLEL!  -----
C
               IF((.NOT.NXT).AND. GOPARR) THEN
                  IPCOUNT= IPCOUNT+1
                  IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 200
               END IF
C
               IANDJ = II.EQ.JJ
               J1 = KSTART(JJ)
               J2 = J1+KNG(JJ)-1
               JPMIN = J1
               JPMAX = J2
               JCNTR = KATOM(JJ)
               JMIN = KMIN(JJ)
               JMAX = KMAX(JJ)
               LOCJ = KLOC(JJ)-JMIN-LOCIJ
               JAMIN = JMIN
               JAMAX = JMAX
               JJMAX = 1
               IF(JMIN.EQ.1 .AND. JMAX.EQ.4) JJMAX = 2
               DO 190 JJJ=1,JJMAX
                  IF(JJMAX.EQ.2) THEN
                     IF(JJJ.EQ.1) THEN
                        JAMIN = 1
                        JAMAX = 1
                     ELSE
                        IF(IANDJ .AND. IAMIN.EQ.1) GO TO 190
                        JAMIN = 2
                        JAMAX = 4
                     END IF
                  END IF
                  DO 150 JG=JPMIN,JPMAX
                     IF(JAMIN.LE.35) T01= CG(JG)
                     IF(JAMIN.LE.20) T01= CF(JG)
                     IF(JAMIN.LE.10) T01= CD(JG)
                     IF(JAMIN.LE. 4) T01= CP(JG)
                     IF(JAMIN.EQ. 1) T01= CS(JG)
                     COEFJ(JG-JPMIN+1)= T01
  150             CONTINUE
                  NPB0= IANG(JAMAX)
                  IJMAX= MAX0((IAMAX-IAMIN+1)*(JAMAX-JAMIN+1),60)
C
                  CALL VCLR(G,1,IJMAX)
C
                  CANDB= ICNTR.EQ.JCNTR
                  R12= ZER
                  DO 160 M=1,3
C CENTER C
                     P12(M,1)= C(M,ICNTR)
C CENTER B
                     P12(M,2)= C(M,JCNTR)
  160             R12= R12+(P12(M,2)-P12(M,1))*(P12(M,2)-P12(M,1))
C
C NOW LOOP OVER EACH CENTER WITH AN ECP POTENTIAL
C
                  NPC= NPC0
                  NPB= NPB0
C N+N' IS THE SUM OF THE ANGULAR MOMENTUM PLUS 1 TO INDEX ARRAYS
                  NPNP= NPC+NPB-1
                  DO 170 IKCNTR=NATST,NATED
                     KCNTR= IKCNTR
C CENTER A
                     ACO(1)= C(1,KCNTR)
                     ACO(2)= C(2,KCNTR)
                     ACO(3)= C(3,KCNTR)
                     IF(LPSKIP(KCNTR).EQ.1) GO TO 170
                     LLMX = LMAX(KCNTR)+1
                     KF1 = KFRST(KCNTR,1)
                     KL1 = KLAST(KCNTR,1)
                     CANDA= ICNTR.EQ.KCNTR
                     AANDB= KCNTR.EQ.JCNTR
C
                     CALL ECPCBA(CANDA,AANDB,ICAB,IPOW)
C
C SET UP TABLES OF THE POWERS OF THE CARTESIAN DISTANCES (CAX, CAY ...)
C FOR LATER USE. PASS IN THE MAXIMUM ANGULAR MOMENTUM FOR I AND J
C USE MAX TO INDEX IANG TO MAKE SURE WE GET THE MAX (L SHELLS)
C
                     IF(ICNTR.NE.KCNTR .OR. KCNTR.NE.JCNTR)
     *               CALL ECPPWR(IPOW,NPC0,NPB0)
C
                     IF(ICAB.EQ.1) THEN
                        CALL ECPAA1(NPNP,FPQR,COEFI,COEFJ,DCOEF4,G)
                     ELSE IF(ICAB.EQ.2 .OR. ICAB.EQ.3) THEN
                        CALL ECPRA2(ICAB,NPNP,FP,COEFI,COEFJ,
     *                              DCOEF1,JFST1,LBECP1,DCOEF4,
     *                              DCOEF2,JFST2,LBECP2,G,
     *                              ZLM,LMF,LMX,LMY,LMZ)
                     ELSE IF(ICAB.EQ.4) THEN
                        IC4C= 0
                        CALL ECPDRA(IC4C,NPNP,FP,FQ,COEFI,COEFQ,COEFJ,
     *                              DCOEF1,JFST1,LBECP1,
     *                              DCOEF2,JFST2,LBECP2,G,
     *                              ZLM,LMF,LMX,LMY,LMZ)
                     END IF
  170             CONTINUE
C
C END OF KCNTR ECP POTENTIAL CENTER LOOP
C
C        -----  STORE ECP INTEGRALS IN ARRAY.                    -----
C
                  MMAX= JAMAX
                  NN= 1
                  DO I=IAMIN,IAMAX
                     IN= IA(LOCI+I)+LOCJ
                     IF(IANDJ) MMAX= I
                     DO J=JAMIN,MMAX
                        HECP(IN+J)= G(NN)
                        NN= NN+1
                     ENDDO
                  ENDDO
C
C END OF BOTH PAIRS OF SHELL LOOPS
C
  190          CONTINUE
  200       CONTINUE
C
  230    CONTINUE
  240 CONTINUE
C
C SYMMETRIZE
C
      IF(INAT.EQ.0) CALL SYMH(HECP,H0,IA)
C
C     -----  ALL DONE, MODIFY BARE NUCLEUS HAMILTONIAN -----
C
      IF(GOPARR) THEN
         IF(NXT) CALL DDI_DLBRESET()
         CALL DDI_GSUMF(900,HECP,L2)
      END IF
C
      IF(DBUG) THEN
         WRITE(IW,*) 'ECP INTEGRAL MODIFICATIONS'
         CALL PRTRI(HECP,NUM)
      END IF
C
      IF(INAT.EQ.0) THEN
        CALL DAREAD(IDAF,IODA,H0,L2,11,0)
        CALL VADD(H0,1,HECP,1,H0,1,L2)
        CALL DAWRIT(IDAF,IODA,H0,L2,11,0)
        IF(ISAVE.EQ.1) CALL DAWRIT(IDAF,IODA,H0,LL2,311,0)
C
        IF(DBUG) THEN
           WRITE(IW,*) '1E INTEGRALS, WITH ECP CHANGES'
           CALL PRTRI(H0,NUM)
        END IF
      END IF
C
      RETURN
      END
C*MODULE ECP     *DECK ECPPWR
      SUBROUTINE ECPPWR(ISWTCH,NI,NJ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C        -----  ROUTINE FINDS THE   NECESSARY POWERS OF THE     -----
C        -----  CARTESIAN COMPONENT VECTORS SEPARATING THE      -----
C        -----  BASIS FUNCTION AND OPERATOR CENTERS.            -----
C THIS ROUTINE CALCULATES (CAX)**(NX-KX) WHERE KX RANGES FROM 0 TO NX
C THUS WE REALLY NEED CAX**(N) WHERE N RANGES FROM 0 TO NX. THESE POWERS
C ARE COMPUTED AND STORED IN THE ACX... ARRAYS WITH ACX(N) BEING
C ACX**(N-1) POWER. NI AND NJ ARE THE MAXIMUM ANGULAR MOMENTUM FOR THE
C I AND J SHELL. HERE WE PRECOMPUTE ALL REQUIRED POWERS ONCE FOR EACH
C COMBINATION OF SHELLS.
C
      COMMON /ECP1  / X01,CAX,CAY,CAZ,CA,XCA,YCA,ZCA,
     *                X02,BAX,BAY,BAZ,BA,XBA,YBA,ZBA,
     *                PHASE,DAX,DAY,DAZ,DA,XDA,YDA,ZDA,XINT,KCNTR
C 7 FOR G HESSIANS
      COMMON /ECP3  / ACX(7),ACY(7),ACZ(7),ABX(7),ABY(7),ABZ(7)
C
      PARAMETER (ONE=1.0D+00)
C FIRST SET THE FIRST ELEMENT TO ONE (0 POWER)
      ACX(1)= ONE
      ACY(1)= ONE
      ACZ(1)= ONE
      IF(ISWTCH.LT.1) THEN
         IF(NI.GT.1) THEN
            ACX(2)=-CAX
            ACY(2)=-CAY
            ACZ(2)=-CAZ
            IF(NI.GT.2) THEN
               DO 110 II=3,NI
                  ACX(II)=-CAX*ACX(II-1)
                  ACY(II)=-CAY*ACY(II-1)
                  ACZ(II)=-CAZ*ACZ(II-1)
  110          CONTINUE
            END IF
         END IF
      END IF
      ABX(1)= ONE
      ABY(1)= ONE
      ABZ(1)= ONE
      IF(ISWTCH.GT.(-1)) THEN
         IF(NJ.GT.1) THEN
            ABX(2)=-BAX
            ABY(2)=-BAY
            ABZ(2)=-BAZ
            IF(NJ.GT.2) THEN
               DO 120 II=3,NJ
                  ABX(II)=-BAX*ABX(II-1)
                  ABY(II)=-BAY*ABY(II-1)
                  ABZ(II)=-BAZ*ABZ(II-1)
  120          CONTINUE
            END IF
         END IF
      END IF
      RETURN
      END
C*MODULE ECP     *DECK ERRF
      DOUBLE PRECISION FUNCTION ERRF(Y)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C        -----  ROUTINE EVALUATES THE ERROR FUNCTION.           -----
C        -----  THE ERROR FUNCTION IS AN ODD-PARITY FUNCTION.   -----
C        -----  EVALUATE   BY A PIECEWISE CHEBYSHEV POLYNOMIAL  -----
C        -----  IF NEEDED USE THE ASYMPTOTIC VALUE.             -----
C
      COMMON /ERRFCM/ C(142),IFRST(20),ILAST(20),H
C
      PARAMETER (ZER=0.0D+00, ONE=1.0D+00)
C
      X= Y
      S= ONE
      IF(X.GE.ZER) GO TO 105
      X=-Y
      S=-ONE
  105 IF(X.GE.4.86D+00) GO TO 115
      XN = X/H
      NX = INT(XN)+1
      I1 = IFRST(NX)
      IL = ILAST(NX)
      T= C(IL)
      DO 110 K=IL-1,I1,-1
  110 T= C(K)+T*X
      ERRF  = T*S
      GO TO 999
  115 ERRF  = S
  999 CONTINUE
      RETURN
      END
C*MODULE ECP     *DECK FA
      DOUBLE PRECISION FUNCTION FA(N,ZETA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C ROUTINE EVALUATES THE TYPE 1 RADIAL INTEGRAL FOR THE SPECIAL
C CASE OF K=-2(ZETA*CA+ZETAB*CB)=0. THUS, Q(N,LAMDA,KR) REDUCES
C TO INTEGRAL FROM ZERO TO INFINITY OF [R**N E**(-ZETA*R**2)]
C GOOD THROUGH FOR N=0 - N=17
C
      DIMENSION GAMMO(9),GAMME(9)
      DATA GAMMO/0.5D+00,0.5D+00,1.0D+00,3.0D+00,12.0D+00,60.0D+00,
     *           360.0D+00,2520.0D+00,20160.0D+00/
      DATA GAMME/0.5D+00,0.25D+00,0.375D+00,0.9375D+00,3.28125D+00,
     *           14.765625D+00,81.2109375D+00,527.87109375D+00,
     *           3959.033203125D+00/
      DATA SQRPI/1.772453850905D+00/
C
      K= N/2
      IF(MOD(N,2).NE.0) THEN
C N = ODD TERMS
         T= GAMMO(K+1)
      ELSE
C N = EVEN TERMS
         T= GAMME(K+1)*SQRPI*SQRT(ZETA)
      END IF
      FA= T/ZETA**(K+1)
      RETURN
      END
C*MODULE ECP     *DECK FIECP
      SUBROUTINE FIECP(FIP,ALFI,NLP,NPNP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION FIP(*)
C
      COMMON /FICMN / ALF,XI,XP0,XP1
C
      DIMENSION FIT(19,2)
C
      PARAMETER (ZP5=0.5D+00)
C
C        -----  ROUTINE PREPARES TABLE OF I- INTEGRALS BY USE   -----
C        -----  OF A RECURSION RELATIONSHIP AND THEN IT PICKS   -----
C        -----  OUT AMONG THESE THE CURRENTLY REQUIRED ONES.    -----
C        -----  IN THIS MANNER ONE CALCULATES ONLY THE UNIQUE   -----
C        -----  MINMAL NUMBER OF I-INTEGRALS.                   -----
C
C        -----  FIND THE LARGEST N VALUE OCCURING IN BOTH THE   -----
C        -----  EVEN AND ODD LATTICES.                          -----
C
      NEFF=NLP+NPNP
      IF(MOD(NEFF,2).NE.0) THEN
         NEMX= NEFF+1
         NOMX= NEFF
      ELSE
         NEMX= NEFF
         NOMX= NEFF+1
      END IF
      IF(NEFF.EQ.0) NOMX=-1
C        -----  PREPARE VARIABLES.                              -----
      YI= XI*ZP5
      A1= XI*ALF
      A2= XI*YI
C        -----  GENERATE A TABLE OF ELEMENTS OF EVEN LATTICE    -----
      IF(NEMX.LT.0) GO TO 115
      FIT(1,1)= FSI0(0)*YI
      IF(NEMX.EQ.0) GO TO 115
      FIT(2,2)= FSI1(1)*YI*YI
      IF(NEMX.EQ.1) GO TO 115
      T1= A2
      DO 110 N=2,NEMX,2
         T2= T1-A2
         FIT(N+1,1)= A1*FIT(N  ,2)+T1*FIT(N-1,1)
         FIT(N+2,2)= A1*FIT(N+1,1)+T2*FIT(N  ,2)
  110 T1= T1+A2+A2
  115 CONTINUE
C        -----  GENERATE TABLE OF ELEMENTS OF ODD-LATTICE       -----
      IF(NOMX.LT.0) GO TO 125
      FIT(1,2)= FSI1(0)*YI
      IF(NOMX.EQ.0) GO TO 125
      FIT(2,1)= FSI0(1)*YI*YI
      IF(NOMX.EQ.1) GO TO 125
      T1= A2+A2
      DO 120 N=2,NOMX,2
         T2= T1-A2-A2-A2
         FIT(N+1,2)= A1*FIT(N  ,1)+T2*FIT(N-1,2)
         FIT(N+2,1)= A1*FIT(N+1,2)+T1*FIT(N  ,1)
  120 T1= T1+A2+A2
  125 CONTINUE
C        -----  RETRIEVE REQUIRED INTEGS FROM TABLES.           -----
      LMAX= 1
      NK=0
      NM1K=0
      DO 150 N=1,NPNP+1
CC       NK=(N*(N-1))/2
         DO 130 L=1,LMAX
  130    FIP(NK+L)= FIT(N+NLP,L)
         LMAX= 2
C APPLY MODIFIED SPHERICAL BESSEL FUNCTION RECURSION RELATION TO
C GET INTEGRALS FOR ALL HIGHER TERMS
         IF(N.LE.2) GO TO 145
CC       NM1K=((N-1)*(N-2))/2
         T01 = ALFI+ALFI+ALFI
         DO 140 L=3,N
            FIP(NK+L)= FIP(NK-2+L)-T01*FIP(NM1K-1+L)
  140    T01 = T01 +ALFI+ALFI
  145    NM1K=NK
  150 NK=NK+N
      RETURN
      END
C*MODULE ECP     *DECK FJFORM
      SUBROUTINE FJFORM(RHO,SGA,SGB,TAU,NMAX)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION RHO(15),SGA(15),SGB(15),TAU(15)
C
C ROUTINE FORMS ALL OF THE NEEDED RADIAL INTEGRALS AND
C STORES THEM IN RHO...TAU SO THAT THE FINAL INTEGRAL MAY
C BE FORMED ELSEWHERE
      COMMON /FJCMN / ALEF,BET,XI,XPLS,XMNS,XP
      COMMON /FJNEW / XKA,XKB,GAM1,GAM2,A1,A2,C
C
      PARAMETER (ZER=0.0D+00, ZP5=0.5D+00, ONE=1.0D+00)
C FACTORIAL(I)*I
      DIMENSION FACTI(13),B1(15),B2(15),C1(15),C2(15)
      DATA FACTI/1.0D+00,4.0D+00,1.8D+01,9.6D+01,6.0D+02,4.32D+03,
     *           3.528D+04,3.2256D+05,3.26592D+06,3.6288D+07,
     *           4.390848D+08,5.7480192D+09,8.09512704D+10/
      DATA SQPI/1.772453850905D+00/
C
      SQPIDC= SQPI*XI
      B1(1) = SQPIDC*GAM1
      C1(1) = SQPIDC*GAM1*ERRF(A1)
      B2(1) = SQPIDC*GAM2
      C2(1) = SQPIDC*GAM2*ERRF(A2)
      TWSQPI= SQPI+SQPI
      B1(2) = TWSQPI*GAM1*DAWERF(A1)
      C1(2) = TWSQPI*GAM1*DAWF(A1)
      B2(2) = TWSQPI*GAM2*DAWERF(A2)
      C2(2) = TWSQPI*GAM2*DAWF(A2)
C
C USE RECURRENCE FORMULAS TO GENERATE ALL OF THE OTHER
C B AND C INTEGRALS. NOTE: GAMMA CANCELS OUT THE EXP IN THE
C RECURRENCE FORMULA TO LEAVE XP. ERRORS MAY ACCUMULATE IF
C GAMMA AND THE EXP ARE USED INSTEAD.
C
      TWC= C+C
      XKAPKB= XKA+XKB
      XKAMKB= XKA-XKB
      XPP= XP
      XPM= XP
      FAP= ZER
      DEN= ONE
      DO 110 I=3,NMAX
         XPP= XPP*XKAPKB
         XPM= XPM*XKAMKB
CC       FAP=(1+((-1)**I))*0.25D+00
CC       FAM=(1-((-1)**I))*0.25D+00
         FAM= ZP5-FAP
         T01= FAP/FACTI(I-2)
         T02= FAM/FACTI(I-2)
         T03= ONE/DEN
         B1(I)= XPP*T01-(TWC*B1(I-2)-XKAPKB*C1(I-1))*T03
         C1(I)= XPP*T02-(TWC*C1(I-2)-XKAPKB*B1(I-1))*T03
         B2(I)= XPM*T01-(TWC*B2(I-2)-XKAMKB*C2(I-1))*T03
         C2(I)= XPM*T02-(TWC*C2(I-2)-XKAMKB*B2(I-1))*T03
      FAP= FAM
  110 DEN= DEN+ONE
C
C STORE THE INTEGRALS IN RHO...
C
      DO 120 I=1,NMAX
         RHO(I)= B1(I)-B2(I)
         SGA(I)= C1(I)+C2(I)
         SGB(I)= C1(I)-C2(I)
         TAU(I)= B1(I)+B2(I)
  120 CONTINUE
      RETURN
      END
C*MODULE ECP     *DECK FJECP
      SUBROUTINE FJECP(FJPQ,ALFI,BETI,NLP,NPNP,LMAX,LEMX)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C        -----  ROUTINE PREPARES TABLE OF J-INTEGS BY USE OF    -----
C        -----  A RECURSION RELATIONSHIP AND THEN IT PICKS UP   -----
C        -----  FROM AMONG THESE THE ONES THAT ARE CURRENTLY    -----
C        -----  NEEDED. THUS
C        -----  ONE ONLY CALCULATES THE UNIQUE NUMBER.          -----
C
      DIMENSION FJPQ(11,11,11)
C
      COMMON /FJCMN / ALEF,BET,XI,XPLS,XMNS,XP
      COMMON /FJNEW / XKA,XKB,GAM1,GAM2,A1,A2,C
C
      DIMENSION RHO(15),SGA(15),SGB(15),TAU(15),FJT(19,2,2)
C
      PARAMETER (ZER=0.0D+00, ZP5=0.5D+00, ONE=1.0D+00, ABLIM=1.0D-01)
      PARAMETER (ZP1=0.1D+00, FIV=5.0D+00, SIX=6.0D+00, TEN=10.0D+00)
C
C        -----  FIND THE LARGEST N VALUES OCCURING IN BOTH THE  -----
C        -----  EVEN AND ODD LATTICES.                          -----
C
      JEFF= NLP+NPNP-1
      NEMX= JEFF
      NOMX= JEFF
      IF(JEFF.EQ.0) NOMX=-1
      LOMX= LEMX
C        -----  IF THE J-INTEGRALS ARE GOING TO BE CALCULATED   -----
C        -----  BY A POWER SERIES, THEN PREPARE A TABLE OF      -----
C        -----  I INTEGS.                                       -----
      IF(ALEF*BET.GT.ABLIM) THEN
         XKAKB = XKA*XKB
         ONEDAB= ONE/XKAKB
         ONDAB2= ONEDAB*ONEDAB
C        IF(2*LEMX+3.GT.15) THEN
         IF(LEMX.GT.6) THEN
            WRITE(6,*) 'LEMX TOO BIG, LEMX= ',LEMX
            CALL ABRT
         END IF
         DO 110 II=1,15
            RHO(II)= ZER
            SGA(II)= ZER
            SGB(II)= ZER
            TAU(II)= ZER
  110    CONTINUE
         CALL FJFORM(RHO,SGA,SGB,TAU,2*LEMX+3)
         FJT(1,1,1)= RHO(3)*ONEDAB
         FJT(2,1,1)= RHO(2)*ONEDAB
         FJT(1,2,1)=(SGB(3)-RHO(4)*ONEDAB*XKB)*ONEDAB
         FJT(2,2,1)=(SGB(2)-RHO(3)*ONEDAB*XKB)*ONEDAB
         FJT(1,1,2)=(SGA(3)-RHO(4)*ONEDAB*XKA)*ONEDAB
         FJT(2,1,2)=(SGA(2)-RHO(3)*ONEDAB*XKA)*ONEDAB
         FJT(1,2,2)=(RHO(5)-XKA*SGB(4)-XKB*SGA(4)+XKAKB*TAU(3))*ONDAB2
         FJT(2,2,2)=(RHO(4)-XKA*SGB(3)-XKB*SGA(3)+XKAKB*TAU(2))*ONDAB2
      ELSE
         CALL SITABL(LEMX,LOMX)
         FJT(1,1,1)= FJPS(0,0,0)
         FJT(2,1,1)= FJPS(1,0,0)
         FJT(1,2,1)= FJPS(0,1,0)
         FJT(2,2,1)= FJPS(1,1,0)
         FJT(1,1,2)= FJPS(0,0,1)
         FJT(2,1,2)= FJPS(1,0,1)
         FJT(1,2,2)= FJPS(0,1,1)
         FJT(2,2,2)= FJPS(1,1,1)
      END IF
      XA= XI*ALEF
      XB= XI*BET
      XX= XI*XI*ZP5
C        -----  GENERATE A TABLE OF ELEMENTS OF THE EVEN        -----
C        -----  LATTICE OF J-INTEGS.                            -----
      IF(NEMX.LT.0) GO TO 125
      IF(NEMX.EQ.0) GO TO 125
      IF(NEMX.EQ.1) GO TO 125
      T2= ZER
      DO 120 N=2,NEMX,2
         T1= T2+XX
         T3= T2-XX-XX-XX
         FJT(N+1,1,1)= XA*FJT(N  ,2,1)+XB*FJT(N  ,1,2)+T1*FJT(N-1,1,1)
         FJT(N+1,2,2)= XA*FJT(N  ,1,2)+XB*FJT(N  ,2,1)+T3*FJT(N-1,2,2)
         FJT(N+2,2,1)= XA*FJT(N+1,1,1)+XB*FJT(N+1,2,2)+T2*FJT(N  ,2,1)
         FJT(N+2,1,2)= XA*FJT(N+1,2,2)+XB*FJT(N+1,1,1)+T2*FJT(N  ,1,2)
  120 T2= T2+XX+XX
  125 CONTINUE
C        -----  GENERATE TABLE OF ELEMENTS OF THE ODD LATTICE   -----
C        -----  J-INTEGS.                                       -----
      IF(NOMX.LT.0) GO TO 135
      IF(NOMX.EQ.0) GO TO 135
      IF(NOMX.EQ.1) GO TO 135
      T2=-XX
      DO 130 N=2,NOMX,2
         T1= T2-XX
         T3= T2+XX+XX+XX
         FJT(N+1,2,1)= XA*FJT(N  ,1,1)+XB*FJT(N  ,2,2)+T2*FJT(N-1,2,1)
         FJT(N+1,1,2)= XA*FJT(N  ,2,2)+XB*FJT(N  ,1,1)+T2*FJT(N-1,1,2)
         FJT(N+2,1,1)= XA*FJT(N+1,2,1)+XB*FJT(N+1,1,2)+T3*FJT(N  ,1,1)
         FJT(N+2,2,2)= XA*FJT(N+1,1,2)+XB*FJT(N+1,2,1)+T1*FJT(N  ,2,2)
  130 T2= T2+XX+XX
  135 CONTINUE
C        -----  RETRIEVE THE REQUIRED J-INTEGS FROM TABLES.     -----
C FIRST TAKE CARE OF KA=0 BLOCK
C
      IF(2.LE.LEMX .AND. LEMX.LE.6) THEN
         IF(ALEF*BET.GT.ABLIM) THEN
            IF(NLP.GT.2) THEN
               WRITE(6,*) 'NLP TOO BIG!, ',NLP,' LEMX = ',LEMX
               CALL ABRT
            END IF
         END IF
      ELSEIF(LEMX.GT.6) THEN
         WRITE(6,*) 'ERROR LEMX = ',LEMX,' IS TOO LARGE IN FJECP'
         CALL ABRT
      END IF
C
      N= NLP
      DO 150 IP= 1,LEMX+1
         IF(IP.LT.3) THEN
            FJPQ(IP,IP,1)= FJT(N+1,IP,IP)
            GO TO 150
         ELSEIF(ALEF*BET.LE.ABLIM) THEN
            FJPQ(IP,IP,1)= FJPS(N,IP-1,IP-1)
            GO TO 150
         END IF
         IF(IP.EQ.3) THEN
            XKA2 = XKA*XKA
            XKB2 = XKB*XKB
            XKAKB2= XKAKB*XKAKB
            DU4= ONDAB2
            QA2= 3.0D+00*XKA2
            QB2= 3.0D+00*XKB2
            T02= QA2+QB2
            T00= 9.0D+00
      DU0 = T00*RHO( 7-N)+T02*RHO( 5-N)+XKAKB2*RHO( 3-N)
      DU1 = T00*SGA( 6-N)+QA2*SGA( 4-N)
      DU2 = T00*SGB( 6-N)+QB2*SGB( 4-N)
      DU3 = T00*TAU( 5-N)
         ELSEIF(IP.EQ.4) THEN
            QA2= 1.5D+01*XKA2
            QB2= 1.5D+01*XKB2
            T02= QA2+QB2
            XRE= SIX*XKAKB2
            T00= 2.25D+02
      DU0 = T00*RHO( 9-N)+(T02*RHO( 7-N)+RHO( 5-N)*XRE)*SIX
      DU1 = T00*SGA( 8-N)+(T02+FIV* QA2)*SGA( 6-N)+SGA( 4-N)*XRE
      DU2 = T00*SGB( 8-N)+(T02+FIV* QB2)*SGB( 6-N)+SGB( 4-N)*XRE
      DU3 = T00*TAU( 7-N)+ T02*TAU( 5-N)+TAU( 3-N)*XKAKB2
         ELSEIF(IP.EQ.5) THEN
            XKA4 = XKA2*XKA2
            XKB4 = XKB2*XKB2
            QA2= 4.5D+01*XKA2
            QB2= 4.5D+01*XKB2
            T02= QA2+QB2
            T03= 1.05D+03*(XKA2+XKB2)
            QA4= 1.05D+02*XKA4
            QB4= 1.05D+02*XKB4
            T04= QA4+QB4
            XRE= TEN*XKAKB2
            T00= 1.1025D+04
      DU0 = T00*RHO(11-N)+1.05D+02*T02*RHO( 9-N)+T04*RHO( 7-N)+
     2     (2.025D+03*RHO( 7-N)+T02*RHO( 5-N)+RHO( 3-N)*XKAKB2)*XKAKB2
      DU1 = T00*SGA(10-N)+(T03+3.675D+03*XKA2)*SGA( 8-N)+
     2      QA4*SGA( 6-N)+(4.5D+01*SGA( 6-N)+XKA2*SGA( 4-N))*XRE
      DU2 = T00*SGB(10-N)+(T03+3.675D+03*XKB2)*SGB( 8-N)+
     2      QB4*SGB( 6-N)+(4.5D+01*SGB( 6-N)+XKB2*SGB( 4-N))*XRE
      DU3 = T00*TAU( 9-N)+T03*TAU( 7-N)+TEN*TAU( 5-N)*XRE
         ELSEIF(IP.EQ.6) THEN
            QA2= 4.2D+02*XKA2
            QB2= 4.2D+02*XKB2
            T02= QA2+QB2
            T03= 9.9225D+04*(XKA2+XKB2)
            QA4= 9.45D+02*XKA4
            QB4= 9.45D+02*XKB4
            T04= QA4+QB4
            XRE= 1.5D+01*XKAKB2
            T00= 8.93025D+05
      DU0 = T00*RHO(13-N)+4.0D+00*T03*RHO(11-N)+1.5D+01*T04*RHO( 9-N)+
     2     (1.176D+04*RHO( 9-N)+T02*RHO( 7-N)+RHO( 5-N)*XRE)*XRE
      DU1 = T00*SGA(12-N)+(T03+7.0875D+02*QA2)*SGA(10-N)+
     2      (T04+1.4D+01*QA4)*SGA( 8-N)+
     3     (4.41D+04*SGA( 8-N)+(T02+2.75D+00*QA2)*SGA( 6-N)+
     4      SGA( 4-N)*XRE)*XKAKB2
      DU2 = T00*SGB(12-N)+(T03+7.0875D+02*QB2)*SGB(10-N)+
     2      (T04+1.4D+01*QB4)*SGB( 8-N)+
     3      (4.41D+04*SGB( 8-N)+(T02+2.75D+00*QB2)*SGB( 6-N)+
     4      SGB( 4-N)*XRE)*XKAKB2
      DU3 = T00*TAU(11-N)+T03*TAU( 9-N)+T04*TAU( 7-N)+
     2     (1.1025D+04*TAU( 7-N)+2.5D-01*T02*TAU( 5-N)+
     3      TAU( 3-N)*XKAKB2)*XKAKB2
         ELSEIF(IP.EQ.7) THEN
            XKA6 = XKA4*XKA2
            XKB6 = XKB4*XKB2
            QA2= 3.5D+00*XKA2
            QB2= 3.5D+00*XKB2
            T02= QA2+QB2
            T03= 6.496875D+01*T02
            QA4= 3.5D+00*XKA4
            QB4= 3.5D+00*XKB4
            T04= 3.75D-01*(QA4+QB4)
            T05= 7.7D+00*T04
            QA6= 4.8125D-02*XKA6
            QB6= 4.8125D-02*XKB6
            T06= QA6+QB6
            XRE= ZP1*XKAKB2/SIX
            T00= 5.00259375D+02
      DU0 = T00*RHO(15-N)+T03*RHO(13-N)+T05*RHO(11-N)+T06*RHO( 9-N)+
     2     (6.2015625D+03*RHO(11-N)+
     3      7.875D+01*T02*RHO( 9-N)+T04*RHO( 7-N)+
     4      (7.35D+02*RHO( 7-N)+T02*RHO( 5-N)+RHO( 3-N)*XRE)*XRE)*XRE
      DU1 = T00*SGA(14-N)+(T03-1.66753125D+02*XKB2)*SGA(12-N)+
     2      (T05-2.59875D+00*QB4)*SGA(10-N)+QA6*SGA( 8-N)+
     3     (1.65375D+04*SGA(10-N)+2.1D+02*(T02-6.25D-01*QB2)*SGA( 8-N)+
     4      QA4*SGA( 6-N)+(7.35D+02*SGA( 6-N)+QA2*SGA( 4-N))*XRE
     5     )*XRE*ZP1
      DU2 = T00*SGB(14-N)+(T03-1.66753125D+02*XKA2)*SGB(12-N)+
     2      (T05-2.59875D+00*QA4)*SGB(10-N)+QB6*SGB( 8-N)+
     3     (1.65375D+04*SGB(10-N)+2.1D+02*(T02-6.25D-01*QA2)*SGB( 8-N)+
     4      QB4*SGB( 6-N)+(7.35D+02*SGB( 6-N)+QB2*SGB( 4-N))*XRE
     5     )*XRE*ZP1
      DU3 = T00*TAU(13-N)+1.7325D+01*T02*TAU(11-N)+ZP1*T05*TAU( 9-N)+
     2     (2.1D+02*TAU( 9-N)+T02*TAU( 7-N)+3.5D+00*TAU( 5-N)*XRE
     3     )*XRE*2.1D+00
      DU4 = DU4*2.16D+05
         END IF
         DU4= DU4*ONEDAB
         FJPQ(IP,IP,1)=(DU0-DU1*XKB-DU2*XKA+DU3*XKAKB)*DU4
  150 CONTINUE
      IF(NPNP.LE.1) GO TO 999
      DO 220 IK=2,NPNP
         DO 160 IP=1,2
            DO 160 IQ=1,2
  160    FJPQ(IQ,IP,IK)= FJT(NLP+IK,IP,IQ)
         IF(LMAX.LE.2) GO TO 220
         T01= ALFI+ALFI+ALFI
         DO 170 IP=3,LMAX
            FJPQ( 1,IP,IK)= FJPQ( 1,IP-2,IK)-T01*FJPQ( 1,IP-1,IK-1)
            FJPQ( 2,IP,IK)= FJPQ( 2,IP-2,IK)-T01*FJPQ( 2,IP-1,IK-1)
  170    T01= T01 +ALFI+ALFI
         T02= BETI+BETI+BETI
         DO 180 IQ=3,LMAX
            FJPQ(IQ, 1,IK)= FJPQ(IQ-2, 1,IK)-T02*FJPQ(IQ-1, 1,IK-1)
            FJPQ(IQ, 2,IK)= FJPQ(IQ-2, 2,IK)-T02*FJPQ(IQ-1, 2,IK-1)
  180    T02= T02 +BETI+BETI
C
         T01= ALFI+ALFI+ALFI
         DO 210 IP=3,LMAX
            T02= BETI+BETI+BETI
            IF(IP.EQ.3) GO TO 195
            DO 190 IQ=3,IP-1
               FJPQ(IQ,IP,IK)= FJPQ(IQ,IP-2,IK)-T01*FJPQ(IQ,IP-1,IK-1)
  190       T02= T02 +BETI+BETI
  195       CONTINUE
            DO 200 IQ=IP,LMAX
               FJPQ(IQ,IP,IK)= FJPQ(IQ-2,IP,IK)-T02*FJPQ(IQ-1,IP,IK-1)
  200       T02= T02 +BETI+BETI
  210    T01= T01 +ALFI+ALFI
  220 CONTINUE
C
  999 CONTINUE
      RETURN
      END
C*MODULE ECP     *DECK FJPS
      DOUBLE PRECISION FUNCTION FJPS(N,LALF,LBET)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /FJCMN / ALF,BET,XI,XPLS,XMNS,XP
      COMMON /FSICMN/ SI(30,7)
C
      PARAMETER(LIMAB=10)
C DOUBLE FACTORIAL
      DIMENSION DFCTRL(7)
      DATA DFCTRL/1.0D+00,3.0D+00,15.0D+00,105.0D+00,945.0D+00,
     *            1.0395D+04,1.35135D+05/
C
C        -----  ROUTINE EVALUATES J-INTEGRALS BY USE OF A       -----
C        -----  POWER SERIES.                                   -----
C        -----  SET UP THE PARAMETERS IN THE POWER SERIES       -----
C        -----  DEPENDING ON THE RELATIVE VALUES OF THE         -----
C        -----  ALFA AND BETA VARIABLES.                        -----
C
      IF((ALF-BET).GT.0.0D+00) THEN
         L0= LBET+1
         X = BET
      ELSE
         L0= LALF+1
         X = ALF
      END IF
      T = X**(L0-1)/DFCTRL(L0)
      X = X*X
      L1= LALF+LBET+2-L0
      L2= L0+N
      L3= L0+L0-1
C        -----  START THE POWER SERIES SUM.                     -----
      SUM = T*SI(L2,L1)
      DO 110 K=1,LIMAB
         T= T*X/(2*K*(K+K+L3))
         SUM = SUM+T*SI(K+K+L2,L1)
  110 CONTINUE
      FJPS= SUM*(0.5D+00*XI)**(N+1)
      RETURN
      END
C*MODULE ECP     *DECK FSI
      DOUBLE PRECISION FUNCTION FSI(N)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C        -----  ROUTINE EVALUES THE SCALED I(N,L) INTEGS        -----
C        -----  FROM ANALYTIC FORMULAS,   OR IF NECESSARY       -----
C        -----  TRANSFERS TO A POWER SERIES APPROACH.           -----
C
      COMMON /FICMN / ALFA,XI,XP0,XP1
C
      SAVE DAWFS,ERRFS
C
      PARAMETER (ZP5=0.5D+00, ONE=1.0D+00, OP5=1.5D+00)
      PARAMETER (TWO=2.0D+00, TP5=2.5D+00, THR=3.0D+00, FIV=5.0D+00)
      PARAMETER (SIX=6.0D+00, SP5=7.5D+00, FIF=15.0D+00)
      PARAMETER (ALIM=0.317D+00)
      DATA SQPI/1.772453850905D+00/
      DATA ERRFS/0.0D+00/
C
C     -------------
      ENTRY FSI0(N)
C     -------------
      IF(ALFA.GT.ALIM) GO TO 010
      FSI0= FSIPS(N,0)
      GO TO 030
  010 NP1 = N+1
      IF(NP1.EQ.2) GO TO 020
      DAWFS=DAWF(ALFA)
      FSI0= SQPI*XP1*DAWFS/ALFA
      GO TO 030
  020 FSI0= ERRFS/ALFA
  030 FSI = FSI0
      RETURN
C
C     -------------
      ENTRY FSI1(N)
C     -------------
      IF(ALFA.GT.ALIM) GO TO 100
      FSI1= FSIPS(N,1)
      GO TO 140
  100 NP1 = N+1
      ALF2= ALFA*ALFA
      TLF2= ALF2+ALF2
      IF(NP1-2) 110,120,130
  110 ERRFS=SQPI*XP1*ERRF(ALFA)
      FSI1=(ZP5*ERRFS-XP0*ALFA)/ALF2
      GO TO 140
  120 FSI1= SQPI*XP1*(ALFA-DAWFS)/ALF2
      GO TO 140
  130 FSI1=(TWO*XP0*ALFA+(TLF2-ONE)*ERRFS)/ALF2
  140 FSI = FSI1
      RETURN
C
C     -------------
      ENTRY FSI2(N)
C     -------------
      IF(ALFA.GT.ALIM) GO TO 200
      FSI2= FSIPS(N,2)
      GO TO 250
  200 NP1 = N+1
      ALF2= ALFA*ALFA
      ALF3= ALFA*ALF2
      TLF2= ALF2+ALF2
      GO TO(210,220,230,240),NP1
  210 FSI2= ZP5*SQPI*XP1*(OP5*ALFA-(ALF2+OP5)*DAWFS)/ALF3
      GO TO 250
  220 FSI2=(THR*XP0*ALFA+(ALF2-OP5)*ERRFS)/ALF3
      GO TO 250
  230 FSI2= SQPI*XP1*(ALFA*(TLF2-THR)+THR*DAWFS)/ALF3
      GO TO 250
  240 FSI2=(TWO*XP0*ALFA*(TLF2-THR)+(TLF2*(TLF2-TWO)+THR)*ERRFS)/ALF3
  250 FSI = FSI2
      RETURN
C
C     -------------
      ENTRY FSI3(N)
C     -------------
      IF(ALFA.GT.ALIM) GO TO 300
      FSI3= FSIPS(N,3)
      GO TO 360
  300 NP1 = N+1
      ALF2= ALFA*ALFA
      ALF3= ALFA*ALF2
      ALF4= ALFA*ALF3
      TLF2= ALF2+ALF2
      GO TO(310,320,330,340,350),NP1
  310 FSI3=(TWO*XP0*ALFA*(TLF2+SP5)/SIX+ZP5*(ALF2-TP5)*ERRFS)/ALF4
      GO TO 360
  320 FSI3= ZP5*SQPI*XP1*(ALFA*(TLF2-SP5)+OP5*(TLF2+FIV)*DAWFS)/ALF4
      GO TO 360
  330 FSI3=(TWO*XP0*ALFA*(ALF2-SP5)+(TLF2*(ALF2-THR)+SP5)*ERRFS)/ALF4
      GO TO 360
  340 FSI3= SQPI*XP1*(ALFA*(TLF2*(TLF2-FIV)+FIF)-FIF*DAWFS)/ALF4
      GO TO 360
  350 FSI3=(TWO*XP0*ALFA*(TLF2*(TLF2-ONE-THR)+FIF)+
     *      (TLF2*(TLF2*(TLF2-THR)+SIX+THR)-FIF)*ERRFS)/ALF4
  360 FSI = FSI3
      RETURN
      END
C*MODULE ECP     *DECK FSIPS
      DOUBLE PRECISION FUNCTION FSIPS(N,L)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C        -----  ROUTINE EVALUATES THE SCALED I(N,L) INTEGS      -----
C        -----  BY A POWER SERIES METHOD.                       -----
C
      COMMON /FICMN / ALFA,XI,XP0,XP1
C
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00)
C
      PARAMETER (LIMA=10)
      DIMENSION FCTRL(10),DFCTRL(10)
      DATA  FCTRL/ONE,ONE,2.0D+00,6.0D+00,24.0D+00,120.0D+00,720.0D+00,
     1            5040.0D+00,40320.0D+00,362880.0D+00/
      DATA DFCTRL/ONE,ONE,3.0D+00,15.0D+00,105.0D+00,945.0D+00,
     1            10395.0D+00,135135.0D+00,2027025.0D+00,34459425.0D+00/
      DATA SQPI/1.772453850905D+00/
C
C        -----  TEST WHETHER (N+L) IS EVEN OR ODD INTEGER       -----
      NL= N+L
      LAMBDA=NL/2
      IF(MOD(NL,2).NE.0) THEN
C        -----  LATTICE OF ODD  (N+L) SCALED I(N,L) INTEGS.     -----
         T  =  FCTRL(LAMBDA+1)/DFCTRL(L+2)
         FAC= TWO**NL
      ELSE
C        -----  LATTICE OF EVEN (N+L) SCALED I(N,L) INTEGS.     -----
         T  = DFCTRL(LAMBDA+1)/DFCTRL(L+2)
         FAC=(TWO**LAMBDA)*SQPI
      ENDIF
      L2= L+L+1
      X = ALFA*ALFA
C        -----  SUM THE POWER SERIES.                           -----
      SUM= T
      DO 110 K=1,LIMA
         T =(T*X*(K+K+NL-1))/(K*(K+K+L2))
  110 SUM= SUM+T
      FSIPS = FAC*XP0*(ALFA**L)*SUM
      RETURN
      END
C*MODULE ECP     *DECK RECUR
      SUBROUTINE RECUR(NMIN,NMAX,LMAX,X)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /FSICMN/ SI(30,7)
C
      PARAMETER (TWO=2.0D+00)
C
C        -----  ROUTINE USES RECURSION RELATIONS TO GENERATE    -----
C        -----  A TABLE OF SCALED I(N,L) INTEGRALS.             -----
C
      TWX= X+X
      T01=(NMIN-TWO)*TWO
C
      IF(LMAX.EQ.0) GO TO 100
      IF(LMAX.LT.0) GO TO 800
      IF(LMAX.GT.6) GO TO 800
      GO TO (100,200,300,400,500,600), LMAX
C
  100 CONTINUE
      DO N=NMIN,NMAX,2
         T02= T01+TWO
         SI(N+1,1)= T02*SI(N-1,1)+TWX*SI(N  ,2)
         SI(N+2,2)= T01*SI(N  ,2)+TWX*SI(N+1,1)
         T01= T02+TWO
      ENDDO
      RETURN
C
  200 CONTINUE
      DO N=NMIN,NMAX,2
         T02= T01+TWO
         SI(N+1,1)= T02*SI(N-1,1)+TWX*SI(N  ,2)
         SI(N+2,2)= T01*SI(N  ,2)+TWX*SI(N+1,1)
         SI(N+3,3)= T01*SI(N+1,3)+TWX*SI(N+2,2)
         T01= T02+TWO
      ENDDO
      RETURN
C
  300 CONTINUE
      DO N=NMIN,NMAX,2
         T02= T01+TWO
         SI(N+1,1)= T02*SI(N-1,1)+TWX*SI(N  ,2)
         SI(N+2,2)= T01*SI(N  ,2)+TWX*SI(N+1,1)
         SI(N+3,3)= T01*SI(N+1,3)+TWX*SI(N+2,2)
         SI(N+4,4)= T01*SI(N+2,4)+TWX*SI(N+3,3)
         T01= T02+TWO
      ENDDO
      RETURN
C            ABOVE THIS POINT IS EQUIVALENT TO ORIGINAL -RECUR- ROUTINE,
C           BEYOND THIS POINT IS BASED ON A HUNCH.
  400 CONTINUE
      DO N=NMIN,NMAX,2
         T02= T01+TWO
         SI(N+1,1)= T02*SI(N-1,1)+TWX*SI(N  ,2)
         SI(N+2,2)= T01*SI(N  ,2)+TWX*SI(N+1,1)
         SI(N+3,3)= T01*SI(N+1,3)+TWX*SI(N+2,2)
         SI(N+4,4)= T01*SI(N+2,4)+TWX*SI(N+3,3)
         SI(N+5,5)= T01*SI(N+3,5)+TWX*SI(N+4,4)
         T01= T02+TWO
      ENDDO
      RETURN
C
  500 CONTINUE
      DO N=NMIN,NMAX,2
         T02= T01+TWO
         SI(N+1,1)= T02*SI(N-1,1)+TWX*SI(N  ,2)
         SI(N+2,2)= T01*SI(N  ,2)+TWX*SI(N+1,1)
         SI(N+3,3)= T01*SI(N+1,3)+TWX*SI(N+2,2)
         SI(N+4,4)= T01*SI(N+2,4)+TWX*SI(N+3,3)
         SI(N+5,5)= T01*SI(N+3,5)+TWX*SI(N+4,4)
         SI(N+6,6)= T01*SI(N+4,6)+TWX*SI(N+5,5)
         T01= T02+TWO
      ENDDO
      RETURN
C
  600 CONTINUE
      DO N=NMIN,NMAX,2
         T02= T01+TWO
         SI(N+1,1)= T02*SI(N-1,1)+TWX*SI(N  ,2)
         SI(N+2,2)= T01*SI(N  ,2)+TWX*SI(N+1,1)
         SI(N+3,3)= T01*SI(N+1,3)+TWX*SI(N+2,2)
         SI(N+4,4)= T01*SI(N+2,4)+TWX*SI(N+3,3)
         SI(N+5,5)= T01*SI(N+3,5)+TWX*SI(N+4,4)
         SI(N+6,6)= T01*SI(N+4,6)+TWX*SI(N+5,5)
         SI(N+7,7)= T01*SI(N+5,7)+TWX*SI(N+6,6)
         T01= T02+TWO
      ENDDO
      RETURN
C
  800 CONTINUE
      WRITE(6,*) 'RECUR: PASSED AN ILLEGAL VALUE, LMAX=',LMAX
      CALL ABRT
      STOP
      END
C*MODULE ECP     *DECK SITABL
      SUBROUTINE SITABL(LEMAX,LOMAX)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C        -----  ROUTINE PREPARES A TABLE OF SCALED I INTEGS     -----
C        -----  FOR USE IN THE EVALUATION OF THE J INTEGS       -----
C        -----  BY A POWER SERIES.                              -----
C
      COMMON /FICMN / ALFI,XI,XP0,XP1
      COMMON /FJCMN / ALFJ,BETJ,XJ,XPLS,XMNS,XP
      COMMON /FSICMN/ SI(30,7)
C
C        -----  PREPARE PARAMETERS OF SCALED INTEGS DEPENDING   -----
C        -----  ON THE SIGN OF THE ALFA AND BETA VARIABLES.     -----
      IF((ALFJ-BETJ).LE.0.0D+00) THEN
         ALFI= BETJ
         XP1 = XMNS
      ELSE
         ALFI= ALFJ
         XP1 = XPLS
      END IF
      X   = ALFI
      XI  = XJ
      XP0 = XP
C        -----  GENERATE LATTICE OF SCALED I(N,L) INTEGS FOR    -----
C        -----  WHICH (N+L) IS EVEN.                            -----
      SI(1,1)= FSI0(0)
      SI(2,2)= FSI1(1)
      IF(LEMAX.LE.1) GO TO 110
      SI(3,3)= FSI2(2)
      IF(LEMAX.EQ.2) GO TO 110
      SI(4,4)= FSI3(3)
      IF(LEMAX.EQ.3) GO TO 110
      SI(5,5)= FSIPS(4,4)
      IF(LEMAX.EQ.4) GO TO 110
      SI(6,6)= FSIPS(5,5)
      IF(LEMAX.EQ.5) GO TO 110
      SI(7,7)= FSIPS(6,6)
  110 CONTINUE
      CALL RECUR(2,22,LEMAX,X)
C        -----  GENERATE LATTICE OF SCALED I(N,L) INTEGS FOR    -----
C        -----  FOR WHICH (N+L) IS ODD.                         -----
      SI(1,2)= FSI1(0)
      SI(2,1)= FSI0(1)
      SI(3,2)= FSI1(2)
      IF(LOMAX.LE.1) GO TO 120
      SI(4,3)= FSI2(3)
      IF(LOMAX.EQ.2) GO TO 120
      SI(5,4)= FSI3(4)
      IF(LOMAX.EQ.3) GO TO 120
      SI(6,5)= FSIPS(4,5)
      IF(LOMAX.EQ.4) GO TO 120
      SI(7,6)= FSIPS(5,6)
      IF(LOMAX.EQ.5) GO TO 120
      SI(8,7)= FSIPS(6,7)
  120 CONTINUE
      CALL RECUR(3,21,LOMAX,X)
      RETURN
      END
C*MODULE ECP     *DECK ZFN
      SUBROUTINE ZFN(ZFNLM,LMAX,ZLM,LMF,LMX,LMY,LMZ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION ZFNLM(125),ZLM(*),LMF(*),LMX(*),LMY(*),LMZ(*)
C
C        -----  THIS ROUTINE EVALUATES THE REAL SPHERICAL
C        -----  HARMONICS GIVEN THE VALUE OF THE THREE
C        -----  CARTESIAN COMPONENTS OF A UNIT VECTOR.
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /ZFNCM / X,Y,Z
C
      PARAMETER (ZER=0.0D+00)
C
      IF(LMAX.LE.10) THEN
         DO 130 L=0,LMAX
            ID= L*(L+1)+L+1
            DO 120 M=-L,L
               IMN= LMF(ID)
               IMX= LMF(ID+1)-1
               SUM= ZER
               DO 110 I=IMN,IMX
                  DUMMY = ZLM(I)
                  IF(LMX(I).GT.0) DUMMY= DUMMY*(X**LMX(I))
                  IF(LMY(I).GT.0) DUMMY= DUMMY*(Y**LMY(I))
                  IF(LMZ(I).GT.0) DUMMY= DUMMY*(Z**LMZ(I))
  110          SUM= SUM+DUMMY
               ZFNLM(ID)= SUM
  120       ID= ID-1
  130    CONTINUE
      ELSE
C
C IF EXECUTION REACHS THIS POINT THEN AN ERROR HAS OCCURED SO ABORT
C SINCE ZFN DOES NOT IMPLEMENT FUNCTIONS HIGHER THAN L=10
C
         WRITE(IW,9010) LMAX
         CALL ABRT
      END IF
      RETURN
 9010 FORMAT(' ERROR: ZFN MAX LAMDA=10, YOU REQUESTED LAMDA = ',I3)
      END
