CHARMM Element source/gamint/blur.src 1.1
##IF GAMESS
      SUBROUTINE BLURIN(H,LL2)
C-----------------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DOUBLE PRECISION NINE,MOROKM
C
      LOGICAL IANDJ,NORM,DOUBLE,GOPARR,DSKWRK,MASWRK,QBARI,QSHORT
C
      DIMENSION H(LL2)
      DIMENSION VBLK(225),FT(225),DIJ(225),zvblk(225),
     *          IJX(225),IJY(225),IJZ(225),
     *          XIN(125),YIN(125),ZIN(125),
     *          IX(35),IY(35),IZ(35),JX(35),JY(35),JZ(35)
C
      INTEGER MXSH,MXGTOT,MXATM,MXCHRM
      PARAMETER (MXSH=1000, MXGTOT=5000, MXATM=500)
      PARAMETER (MXCHRM=25120)
C
      COMMON /CHMGMS/ XCHM(MXCHRM),YCHM(MXCHRM),ZCHM(MXCHRM),
     *                DXELMM(MXCHRM),DYELMM(MXCHRM),DZELMM(MXCHRM),
     *                QCHM(MXCHRM),NCHMAT,KCHRMM,KGUES
      COMMON /CGAMES/ NGAMES,NQQCHG,QQCHG(MXATM),FQQCHG(MXATM)
      LOGICAL QGMREM,QGMEXG,QINIGM,QBLUCH,QNOGU
      COMMON /GAMESL/QGMREM,QGMEXG,QINIGM,QBLUCH,QNOGU
      INTEGER MAXBLU
      PARAMETER (MAXBLU=25120)
      INTEGER NBLUCH,IBLUCH
      REAL*8 EBLUCH,CBLUCH
      COMMON /BLUR/EBLUCH(MAXBLU),CBLUCH(MAXBLU),IBLUCH(MAXBLU),NBLUCH
      REAL*8 CGBLCH,SGBLCH
      COMMON /MBLUR/ CGBLCH(MAXBLU), SGBLCH(MAXBLU)
C
      COMMON /CBARI/ QX,QY,QZ,CKN,CLN,AK
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NNP,NE,NA,NB,ZAN(MXATM),C(3,MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(400)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ROOT  / XX,U(9),W(9),NROOTS
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS
      COMMON /SCINP / VLAMB,SCREEN
      COMMON /STV   / XINT,YINT,ZINT,TAA,X0,Y0,Z0,
     *                XI,YI,ZI,XJ,YJ,ZJ,NI,NJ
      COMMON /SYMIND/ TOL,II,JJ,LIT,LJT,MINI,MINJ,MAXI,MAXJ,IANDJ
C
C     These are taken from int2a.src:
C-----------------------------------
C                                 NORG needed in SHELLS !
      COMMON /GOUT  / GPOPLE(768),NORG
      COMMON /SHLNOS/ QQ4,LITQ,LJTQ,LKTQ,LLT,LOCIQ,LOCJQ,LOCKQ,LOCLQ,
     *                MINIQ,MINJQ,MINKQ,MINLQ,MAXIQ,MAXJQ,MAXKQ,MAXLQ,
     *                NIJQ,IJQ,KLQ,IJKL
      INTEGER MXGSH, MXG2
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      DIMENSION DDIJ(16*MXG2)
C
C-----------------------------------
      PARAMETER (ZERO=0.0D+00, PT5=0.5D+00, ONE=1.0D+00, TWO=2.0D+00,
     *           THREE=3.0D+00, FIVE=5.0D+00, SEVEN=7.0D+00,
     *           NINE=9.0D+00, ELEVEN=11.0D+00, E999=999.0D+00,
     *           PI212=1.1283791670955D+00, SQRT3=1.73205080756888D+00,
     *           SQRT5=2.23606797749979D+00, SQRT7=2.64575131106459D+00,
     *           RLN10=2.30258D+00,PI=3.14159265358979323844D+00,
     *           PI32=5.56832799683170D+00)
C
      DATA JX / 0, 1, 0, 0, 2, 0, 0, 1, 1, 0,
     *          3, 0, 0, 2, 2, 1, 0, 1, 0, 1,
     *          4, 0, 0, 3, 3, 1, 0, 1, 0, 2,
     *          2, 0, 2, 1, 1/
      DATA IX / 1, 6, 1, 1,11, 1, 1, 6, 6, 1,
     *         16, 1, 1,11,11, 6, 1, 6, 1, 6,
     *         21, 1, 1,16,16, 6, 1, 6, 1,11,
     *         11, 1,11, 6, 6/
      DATA JY / 0, 0, 1, 0, 0, 2, 0, 1, 0, 1,
     *          0, 3, 0, 1, 0, 2, 2, 0, 1, 1,
     *          0, 4, 0, 1, 0, 3, 3, 0, 1, 2,
     *          0, 2, 1, 2, 1/
      DATA IY / 1, 1, 6, 1, 1,11, 1, 6, 1, 6,
     *          1,16, 1, 6, 1,11,11, 1, 6, 6,
     *          1,21, 1, 6, 1,16,16, 1, 6,11,
     *          1,11, 6,11, 6/
      DATA JZ / 0, 0, 0, 1, 0, 0, 2, 0, 1, 1,
     *          0, 0, 3, 0, 1, 0, 1, 2, 2, 1,
     *          0, 0, 4, 0, 1, 0, 1, 3, 3, 0,
     *          2, 2, 1, 1, 2/
      DATA IZ / 1, 1, 1, 6, 1, 1,11, 1, 6, 6,
     *          1, 1,16, 1, 6, 1, 6,11,11, 6,
     *          1, 1,21, 1, 6, 1, 6,16,16, 1,
     *         11,11, 6, 6,11/
      DATA MOROKM/8HMOROKUMA/
C
C     ----- COMPUTE ADDITIONAL H INTEGRALS DUE TO CHARMM FIELD -----
C
      TOL = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
C
C     ----- RESET SOME PARAMETERS FOR MOROKUMA DECOMPOSITIONS -----
C     ISAVE .EQ. 0 : SAVE S, H, AND T TO DAF 12, 11, AND 13
C     ISAVE .EQ. 1 : SAVE S, H, AND T TO DAF 12, 11, AND 13
C                    AND SAVE S AND H TO DAF 312 AND 311
C     NOTE THAT LL2 IS ALWAYS (NUM*NUM+NUM)/2,
C     L1,L2 MAY BE SMALLER THAN USUAL FOR A MONOMER IN A MOROKUMA RUN
C
      IF (RUNTYP.EQ.MOROKM) THEN
         CALL STINT1(ISTART,IEND,JSTART,LOCIJ,NATST,NATED,ISAVE,L1,L2)
      ELSE
         ISTART = 1
         IEND   = NSHELL
         JSTART = 1
         LOCIJ  = 0
         NATST  = NAT+1
         NATED  = NAT+NCHMAT
         ISAVE  = 0
         L1 = NUM
         L2 = (NUM*(NUM+1))/2
      END IF
C
      NINT  = 0
      NSCHWZ= 0
      DENMAX = ZERO
C
C     ----- INTIALIZE PARALLEL -----
C
      IPCOUNT = ME - 1
C
C     This is the loop for external charges. Because
C     the ``blur'' method needs data before the inner loop
C     the loop over the charge
C     centers has to be taken out from the inner loop.
C
C     -NCHMAT- IS NONZERO IF THERE ARE EXTERNAL CHARGES WHICH
C     PERTURB THE SYSTEM, SUCH AS IF CHARMM IS IN USE.  NOTE
C     THAT THERE IS ALSO A NUCLEAR REPULSION TERM WHICH IS NOT
C     INCLUDED HERE, IT IS IN THE CHARMM INTERFACE CODE.
C
C
      IBLUR=1
      DO 460 IC = NATST,NATED
         ZNUC = -QCHM(IC-NAT)
         CX = XCHM(IC-NAT)
         CY = YCHM(IC-NAT)
         CZ = ZCHM(IC-NAT)
C
         QBARI=.FALSE.
         QSHORT=.FALSE.
         AK=ZERO
         IF(QBLUCH.AND.(IBLUCH(IBLUR).EQ.(IC-NAT))) THEN
            AK=PT5/SGBLCH(IBLUR)/SGBLCH(IBLUR)
            ZNUC=-CGBLCH(IBLUR)
            QBARI=.TRUE.
            QSHORT=.FALSE.
CC            QSHORT=.FALSE.
CCC            write(*,*)'BLURIN>iblur,znuc,sgblch,ak,qbari,qshort',
CCC     $           iblur,znuc,sgblch(iblur),ak,qbari,qshort
            IBLUR=IBLUR+1
         ENDIF
C        This is the normalization factor
         CKN=TWO*AK/PI
         CKN=CKN*CKN*CKN
         CKN=SQRT(CKN)
         CKN=SQRT(CKN)
         CLN=ZNUC*CKN
         QX=CX
         QY=CY
         QZ=CZ
CCC         write(*,*)'BLURIN>ak,znuc,cx,cy,cz=',ak,znuc,cx,cy,cz
CCC         write(*,*)'BLURIN>ckn,cln,qx,qy,qz=',ckn,cln,qx,qy,qz
C
C     ----- I SHELL -----
C
CCC         write(*,*)'BLURIN>nshell=',nshell
         DO 720 II = ISTART,IEND
            I = KATOM(II)
            XI = C(1,I)
            YI = C(2,I)
            ZI = C(3,I)
            I1 = KSTART(II)
            I2 = I1+KNG(II)-1
            LIT = KTYPE(II)
            MINI = KMIN(II)
            MAXI = KMAX(II)
            LOCI = KLOC(II)-MINI-LOCIJ
CCC         write(*,*)'BLURIN>xi,yi,zi=',xi,yi,zi
CCC         write(*,*)'BLURIN>ii,i,i1,i2,lit,mini,maxi,loci=',
CCC     $        ii,i,i1,i2,lit,mini,maxi,loci
C      
C     ----- J SHELL -----
C 
            DO 700 JJ = JSTART,II
C
C     ----- GO PARALLEL! -----
C
               IF (GOPARR) THEN
                  IPCOUNT = IPCOUNT + 1
                  IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 700
               END IF
               J = KATOM(JJ)
               XJ = C(1,J)
               YJ = C(2,J)
               ZJ = C(3,J)
               J1 = KSTART(JJ)
               J2 = J1+KNG(JJ)-1
               LJT = KTYPE(JJ)
               MINJ = KMIN(JJ)
               MAXJ = KMAX(JJ)
               LOCJ = KLOC(JJ)-MINJ-LOCIJ
               NROOTS = (LIT+LJT-2)/2+1
               RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
               IANDJ = II .EQ. JJ
C
C     ----- PREPARE INDICES FOR PAIRS OF (I,J) FUNCTIONS
C
               IJ = 0
               MAX = MAXJ
CCC            write(*,*)'BLURIN>xj,yj,zj=',xj,yj,zj
CCC           write(*,*)'BLURIN>jj,j,j1,j2,ljt,minj,maxj,locj=',
CCC     $           jj,j,j1,j2,ljt,minj,maxj,locj
CCC            write(*,*)'BLURIN>nroots,rr,iandj=',nroots,rr,iandj
               DO 160 I = MINI,MAXI
                  NX = IX(I)
                  NY = IY(I)
                  NZ = IZ(I)
                  IF (IANDJ) MAX = I
                  DO 140 J = MINJ,MAX
                     IJ = IJ+1
                     IJX(IJ) = NX+JX(J)
                     IJY(IJ) = NY+JY(J)
                     IJZ(IJ) = NZ+JZ(J)
                     IF (J.LE.1) FT(IJ) = THREE
                     IF ((J.GT.1).AND.(J.LE.4)) FT(IJ) = FIVE
                     IF ((J.GT.4).AND.(J.LE.10)) FT(IJ) = SEVEN
                     IF ((J.GT.10).AND.(J.LE.20)) FT(IJ) = NINE
                     IF (J.GT.20) FT(IJ) = ELEVEN
 140              CONTINUE
 160           CONTINUE
C
               DO 180 I = 1,IJ
                  zVBLK(I) = ZERO
                  VBLK(I) = ZERO
 180           CONTINUE
C
               IF (QBARI.AND.(.NOT.QSHORT)) THEN
C
C     ----- K,L SHELL ----- For BARI ignore the loops and set 
C                           the indeces to be always 1
C
                  IEXCH = 1
                  ISH = II
                  JSH = JJ
                  KSH = 1
                  LSH = 1
                  QQ4 = ONE
C
C        ----- COMPUTE TWO-ELECTRON INTEGRALS ----
C
C     USE HONDO RYS POLYNOMIAL CODE FOR OTHER BLOCKS
C
C        ----- GET INFORMATION ABOUT ISH AND JSH -----
C        ----- FORM PAIRS OF PRIMITIVES FROM ISH AND JSH -----
C        ----- GET INFORMATION ABOUT KSH AND LSH -----
C
                  NORG=0
CCC                  write(*,*)'BLURIN-b-shells>xi,yi,zi=',xi,yi,zi
CCC                  write(*,*)'BLURIN-b-shells>xj,yj,zj=',xj,yj,zj
                  CALL SHELLS(1,ISH,JSH,KSH,LSH)
CCC                  write(*,*)'BLURIN-b-ijprim>xi,yi,zi=',xi,yi,zi
CCC                  write(*,*)'BLURIN-b-ijprim>xj,yj,zj=',xj,yj,zj
                  CALL IJPRIM1(DDIJ)
C     Info which is needed, but is missing because call to
C     SHELLS(2,ISH,JSH,KSH,LSH) would be meaningless for BARI
C     Provide here:
                  IJKL=IJ
C
C        ----- DO INTEGRAL BATCH, SSSS IS A SPECIAL CASE -----
C                                 but general can handle them too!
C
CCC                  IF(IJKL.EQ.100) THEN
CCC                     CALL BARI00(VBLK,DDIJ)
CCC                     write(*,'(a,2f20.10)')'BLURIN>1:ddij,vblk',
CCC     $                    ddij(1),vblk(1)
CCC                  ELSE
CCC                  write(*,*)'BLURIN-b-barign>xi,yi,zi=',xi,yi,zi
CCC                  write(*,*)'BLURIN-b-barign>xj,yj,zj=',xj,yj,zj
                  CALL BARIGN(VBLK,DDIJ)
CCC                  write(*,*)'BLURIN-a-barign>xi,yi,zi=',xi,yi,zi
CCC                  write(*,*)'BLURIN-a-barign>xj,yj,zj=',xj,yj,zj
CCC                     write(*,*)'BLURIN>ij,ddij,vblk=',ij
CCC                     write(*,'(10f12.8)')
CCC     $                    (ddij(i),vblk(i),i=1,ij)
CCC                  END IF
                  GOTO 900
               ENDIF
C
C     ----- I PRIMITIVE
C
               JGMAX = J2
CCC            write(*,*)'BLURIN>ij,ijx,ft=',ij,(ijx(i),i=1,ij),
CCC     $           (ft(i),i=1,ij)
               DO 520 IG = I1,I2
                  AI = EX(IG)
                  ARRI = AI*RR
                  AXI = AI*XI
                  AYI = AI*YI
                  AZI = AI*ZI
                  CSI = CS(IG)
                  CPI = CP(IG)
                  CDI = CD(IG)
                  CFI = CF(IG)
                  CGI = CG(IG)
CCC               write(*,*)'BLURIN>ai,arri,axi,ayi,azi=',
CCC     $              ai,arri,axi,ayi,azi
CCC               write(*,*)'BLURIN>csi,cpi,cdi,cgi=',csi,cpi,cdi,cgi
C
C     ----- J PRIMITIVE
C
                  IF (IANDJ) JGMAX = IG
                  DO 500 JG = J1,JGMAX
                     AJ = EX(JG)
                     AA = AI+AJ
                     AA1 = ONE/AA
                     DUM = AJ*ARRI*AA1
CCC                     write(*,*)'BLURIN>aj,aa,aa1,dum,tol=',
CCC     $                    aj,aa,aa1,dum,tol
                     IF (QSHORT.AND.QBARI.AND.(AK.EQ.ZERO)) GO TO 500
                     IF (DUM .GT. TOL) GO TO 500
                     FAC = EXP(-DUM)
                     CSJ = CS(JG)
                     CPJ = CP(JG)
                     CDJ = CD(JG)
                     CFJ = CF(JG)
                     CGJ = CG(JG)
                     AX = (AXI+AJ*XJ)*AA1
                     AY = (AYI+AJ*YJ)*AA1
                     AZ = (AZI+AJ*ZJ)*AA1
CCC                  write(*,*)'BLURIN>fac,ax,ay,az=',fac,ax,ay,az
CCC                  write(*,*)'BLURIN>csj,cpj,cdj,cgj=',csj,cpj,cdj,cgj
C
C     ----- DENSITY FACTOR
C
CCC                  write(*,*)'BLURIN>mini,maxi,minj,maxj=',
CCC     $                 mini,maxi,minj,maxj
CCC                   write(*,*)'BLURIN>csi,csj,cpi,cpj=',csi,csj,cpi,cpj
                     DOUBLE=IANDJ.AND.IG.NE.JG
                     MAX = MAXJ
                     NN = 0
                     DUM1 = ZERO
                     DUM2 = ZERO
                     DO 220 I = MINI,MAXI
                        IF (I.EQ.1) DUM1=CSI*FAC
                        IF (I.EQ.2) DUM1=CPI*FAC
                        IF (I.EQ.5) DUM1=CDI*FAC
                        IF ((I.EQ. 8).AND.NORM) DUM1=DUM1*SQRT3
                        IF (I.EQ.11) DUM1=CFI*FAC
                        IF ((I.EQ.14).AND.NORM) DUM1=DUM1*SQRT5
                        IF ((I.EQ.20).AND.NORM) DUM1=DUM1*SQRT3
                        IF (I.EQ.21) DUM1=CGI*FAC
                        IF ((I.EQ.24).AND.NORM) DUM1=DUM1*SQRT7
                        IF ((I.EQ.30).AND.NORM) DUM1=DUM1*SQRT5/SQRT3
                        IF ((I.EQ.33).AND.NORM) DUM1=DUM1*SQRT3
                        IF (IANDJ) MAX = I
                        DO 200 J = MINJ,MAX
                           IF (J.EQ.1) THEN
                              DUM2=DUM1*CSJ
                              IF (DOUBLE) THEN
                                 IF (I.LE.1) THEN
                                    DUM2=DUM2+DUM2
                                 ELSE
                                    DUM2=DUM2+CSI*CPJ*FAC
                                 END IF
                              END IF
                           ELSE IF (J.EQ.2) THEN
                              DUM2=DUM1*CPJ
                              IF (DOUBLE) DUM2=DUM2+DUM2
                           ELSE IF (J.EQ.5) THEN
                              DUM2=DUM1*CDJ
                              IF (DOUBLE) DUM2=DUM2+DUM2
                           ELSE IF ((J.EQ.8).AND.NORM) THEN
                              DUM2=DUM2*SQRT3
                           ELSE IF (J.EQ.11) THEN
                              DUM2=DUM1*CFJ
                              IF (DOUBLE) DUM2=DUM2+DUM2
                           ELSE IF ((J.EQ.14).AND.NORM) THEN
                              DUM2=DUM2*SQRT5
                           ELSE IF ((J.EQ.20).AND.NORM) THEN
                              DUM2=DUM2*SQRT3
                           ELSE IF (J.EQ.21) THEN
                              DUM2=DUM1*CGJ
                              IF (DOUBLE) DUM2=DUM2+DUM2
                           ELSE IF ((J.EQ.24).AND.NORM) THEN
                              DUM2=DUM2*SQRT7
                           ELSE IF ((J.EQ.30).AND.NORM) THEN
                              DUM2=DUM2*SQRT5/SQRT3
                           ELSE IF ((J.EQ.33).AND.NORM) THEN
                              DUM2=DUM2*SQRT3
                           END IF
                           NN = NN+1
                           DIJ(NN) = DUM2
 200                    CONTINUE
 220                 CONTINUE
CCC                     write(*,*)'BLURIN>nn,dij=',nn,(dij(i),i=1,nn)
C
C     ----- NUCLEAR ATTRACTION
C
C                    PI212 = TWO/SQRT(PI), PI32 = PI**(3/2)
C
                     DUM = PI212*AA1
                     IF(QBARI.AND.QSHORT)
     $                  DUM=DUM*PI32*CKN*CKN/(TWO*AK*SQRT(AA+TWO*AK))
CCC                     write(*,*)'BLURIN>dum=',dum
                     DO 400 I = 1,IJ
                        DIJ(I) = DIJ(I)*DUM
 400                 CONTINUE
C
                     AAX = AA*AX
                     AAY = AA*AY
                     AAZ = AA*AZ
C
                     IF (QBARI.AND.QSHORT) THEN
                        RHO=AA*TWO*AK/(AA+AK+AK)
                        XX = RHO*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
                     ELSE
                        XX = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
                     ENDIF
CCC                     write(*,*)'BLURIN>xi,yi,zi,xj,yj,zj,cx,cy,cz=',
CCC     $                    xi,yi,zi,xj,yj,zj,cx,cy,cz
CCC                     write(*,*)'BLURIN>aa,rho,xx=',aa,rho,xx
                     IF (NROOTS.LE.3) CALL RT123
                     IF (NROOTS.EQ.4) CALL ROOT4
                     IF (NROOTS.EQ.5) CALL ROOT5
                     MM = 0
                     DO 430 K = 1,NROOTS
CCC                        write(*,*)'BLURIN>root,u(),w()=',k,u(k),w(k)
                        IF(QBARI.AND.QSHORT) THEN
                           UU = RHO*U(K)
                        ELSE
                           UU = AA*U(K)
                        ENDIF
                        WW = W(K)*ZNUC
                        TT = ONE/(AA+UU)
                        TAA = SQRT(TT)
                        X0 = (AAX+UU*CX)*TT
                        Y0 = (AAY+UU*CY)*TT
                        Z0 = (AAZ+UU*CZ)*TT
                        IN = -5+MM
                        DO 420 I = 1,LIT
                           IN = IN+5
                           NI = I
                           DO 410 J = 1,LJT
                              JN = IN+J
                              NJ = J
                              CALL STVINT
CCC                              write(*,*)'BLURIN>xint,yint,zint=',
CCC     $                             xint,yint,zint
                              XIN(JN) = XINT
                              YIN(JN) = YINT
                              ZIN(JN) = ZINT*WW
CCC                              write(*,*)'BLURIN>uu,ww,tt,taa=',
CCC     $                             uu,ww,tt,taa
CCC                              write(*,*)'BLURIN>ni,nj,x0,y0,z0=',
CCC     $                             ni,nj,x0,y0,z0
CCC                              write(*,*)'BLURIN>jn,xin,yin,zin=',
CCC     $                             jn,xin(jn),yin(jn),zin(jn)
 410                       CONTINUE
 420                    CONTINUE
                        MM = MM+25
 430                 CONTINUE
                     DO 450 I = 1,IJ
                        NX = IJX(I)
                        NY = IJY(I)
                        NZ = IJZ(I)
                        DUM = ZERO
                        MM = 0
                        DO 440 K = 1,NROOTS
                           DUM = DUM+XIN(NX+MM)*YIN(NY+MM)*ZIN(NZ+MM)
CCC                           write(*,*)'BLURIN>nx,ny,nz,mm=',nx,ny,nz,mm
CCC                           write(*,*)'BLURIN>k,dum=',k,dum
                           MM = MM+25
 440                    CONTINUE
                        VBLK(I) = VBLK(I) + DUM*DIJ(I)
CCC                        write(*,*)'BLURIN>i,dij,vblk=',i,dij(i),vblk(i)
 450                 CONTINUE
C
C     ----- END OF PRIMITIVE LOOPS -----
C
 500              CONTINUE
 520           CONTINUE
 900           CONTINUE
C
C     ----- COPY BLOCK INTO H-CORE, OVERLAP, AND KINETIC ENERGY MATRICES
C
               MAX = MAXJ
               NN = 0
               DO 620 I = MINI,MAXI
                  LI = LOCI+I
                  IN = (LI*(LI-1))/2
                  IF (IANDJ) MAX = I
                  DO 600 J = MINJ,MAX
                     LJ = LOCJ+J
                     JN = LJ+IN
                     NN = NN+1
                     H(JN) = H(JN) + VBLK(NN)
 600              CONTINUE
 620           CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
 700        CONTINUE
 720     CONTINUE
 460  CONTINUE
C
      RETURN
      END
C
      SUBROUTINE BARI00(GHONDO,DDIJ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C
      LOGICAL IANDJ,KANDL,SAME,OUT
C
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      DIMENSION GHONDO(*),DDIJ(16*MXG2)
C
C     For blur
      COMMON /CBARI/ QX,QY,QZ,CKN,CLN,AK
C
      COMMON /IJGNRL/ A(MXG2),R(MXG2),X1(MXG2),Y1(MXG2),Z1(MXG2),
     *                IJD(225)
      COMMON /GOUT  / GPOPLE(768),NORG
      COMMON /MISC  / IANDJ,KANDL,SAME
      COMMON /SHLINF/ AG(MXGSH),CSA(MXGSH),CPA(MXGSH),CDA(MXGSH),
     *                CFA(MXGSH),CGA(MXGSH),
     *                BG(MXGSH),CSB(MXGSH),CPB(MXGSH),CDB(MXGSH),
     *                CFB(MXGSH),CGB(MXGSH),
     *                CG(MXGSH),CSC(MXGSH),CPC(MXGSH),CDC(MXGSH),
     *                CFC(MXGSH),CGC(MXGSH),
     *                DG(MXGSH),CSD(MXGSH),CPD(MXGSH),CDD(MXGSH),
     *                CFD(MXGSH),CGD(MXGSH),
     *                XI,YI,ZI,XJ,YJ,ZJ,RRI,XK,YK,ZK,XL,YL,ZL,RRK,
     *                NGA,NGB,NGC,NGD
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJ,IJ,KL,IJKL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
C
      PARAMETER (PI252=34.986836655250D+00, PIE4=7.85398163397448D-01,
     *           ZERO=0.0D+00, ONE=1.0D+00)
C
C     SPECIAL SSSS INTEGRAL ROUTINE WHEN USING HONDO INTEGRALS
C
      GGOUT = ZERO
      BK = AK
      RRK=ZERO
      BRRK = ZERO
      BXK = BK*QX
      BYK = BK*QY
      BZK = BK*QZ
      CSK = CKN
      BL = AK
      BB = BK+BL
      BBINV = ONE/BB
      DUM = BL*BRRK*BBINV
      IF (DUM .GT. TOL) GO TO 280
      BBRRK = DUM
      D2 = CLN*CSK*BBINV
      BBX = QX
      BBY = QY
      BBZ = QZ
      SUM = ZERO
      NN = 1
CCC      write(*,*)'BARI00>qx,qy,qz,ckn,cln',qx,qy,qz,ckn,cln
CCC      write(*,*)'BARI00>nij=',nij
      DO 260 N = 1,NIJ
      DUM = BBRRK+R(N)
      IF (DUM .GT. TOL) GO TO 260
      EXPE = EXP(-DUM)
      AA = A(N)
      AB = AA+BB
      DUM = X1(N)-BBX
      XX = DUM*DUM
      DUM = Y1(N)-BBY
      XX = DUM*DUM+XX
      DUM = Z1(N)-BBZ
      XX = DUM*DUM+XX
      X = XX*AA*BB/AB
CCC      write(*,*)'BARI00>aa,ab,xx,bb=',aa,ab,xx,bb
CCC      write(*,*)'BARI00>x=',X
C
      IF (X .GT. 5.0D+00) GO TO 160
      IF (X .GT. 1.0D+00) GO TO 120
      IF (X .GT. 3.0D-07) GO TO 100
      WW1 = 1.0D+00-X/3.0D+00
      GO TO 240
C
  100 CONTINUE
      F1 = ((((((((-8.36313918003957D-08*X+1.21222603512827D-06 )*X-
     +     1.15662609053481D-05 )*X+9.25197374512647D-05 )*X-
     +     6.40994113129432D-04 )*X+3.78787044215009D-03 )*X-
     +     1.85185172458485D-02 )*X+7.14285713298222D-02 )*X-
     +     1.99999999997023D-01 )*X+3.33333333333318D-01
      WW1 = (X+X)*F1+EXP(-X)
      GO TO 240
C
  120 CONTINUE
      IF (X .GT. 3.0D+00) GO TO 140
      Y = X-2.0D+00
      F1 = ((((((((((-1.61702782425558D-10*Y+1.96215250865776D-09 )*Y-
     +     2.14234468198419D-08 )*Y+2.17216556336318D-07 )*Y-
     +     1.98850171329371D-06 )*Y+1.62429321438911D-05 )*Y-
     +     1.16740298039895D-04 )*Y+7.24888732052332D-04 )*Y-
     +     3.79490003707156D-03 )*Y+1.61723488664661D-02 )*Y-
     +     5.29428148329736D-02 )*Y+1.15702180856167D-01
      WW1 = (X+X)*F1+EXP(-X)
      GO TO 240
C
  140 CONTINUE
      Y = X-4.0D+00
      F1 = ((((((((((-2.62453564772299D-11*Y+3.24031041623823D-10 )*Y-
     +     3.614965656163D-09)*Y+3.760256799971D-08)*Y-
     +     3.553558319675D-07)*Y+3.022556449731D-06)*Y-
     +     2.290098979647D-05)*Y+1.526537461148D-04)*Y-
     +     8.81947375894379D-04 )*Y+4.33207949514611D-03 )*Y-
     +     1.75257821619926D-02 )*Y+5.28406320615584D-02
      WW1 = (X+X)*F1+EXP(-X)
      GO TO 240
C
  160 CONTINUE
      IF (X .GT. 15.0D+00) GO TO 200
      E = EXP(-X)
      IF (X .GT. 10.0D+00) GO TO 180
      XINV = ONE/X
      WW1 = (((((( 4.6897511375022D-01*XINV-6.9955602298985D-01)*XINV +
     +     5.3689283271887D-01)*XINV-3.2883030418398D-01)*XINV +
     +     2.4645596956002D-01)*XINV-4.9984072848436D-01)*XINV -
     +     3.1501078774085D-06)*E + SQRT(PIE4*XINV)
      GO TO 240
C
  180 CONTINUE
      XINV = ONE/X
      WW1 = (((-1.8784686463512D-01*XINV+2.2991849164985D-01)*XINV
     +         -4.9893752514047D-01)*XINV-2.1916512131607D-05)*E
     +         + SQRT(PIE4*XINV)
      GO TO 240
C
  200 CONTINUE
      IF (X .GT. 33.0D+00) GO TO 220
      XINV = ONE/X
      E = EXP(-X)
      WW1 = (( 1.9623264149430D-01*XINV-4.9695241464490D-01)*XINV -
     +     6.0156581186481D-05)*E + SQRT(PIE4*XINV)
      GO TO 240
C
  220 WW1 = SQRT(PIE4/X)
C
  240 SUM = SUM+DDIJ(NN)*WW1*EXPE/SQRT(AB)
  260 NN = NN+16
      GGOUT = GGOUT+D2*SUM
 280  CONTINUE
      GHONDO(1) = GGOUT*PI252*QQ4
      RETURN
      END
C
      SUBROUTINE BARIGN(GHONDO,DDIJ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C
      DIMENSION GHONDO(*),DDIJ(*)
C
      LOGICAL IANDJ,KANDL,SAME,OUT,NORM,DOUBLE
C
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      COMMON /CBARI/ QX,QY,QZ,CKN,CLN,AK
C
      COMMON /DENS  / DKL(225),DIJ(225)
      COMMON /IJGNRL/ AA(MXG2),R(MXG2),X1(MXG2),Y1(MXG2),Z1(MXG2),
     *                IJD(225)
      COMMON /MISC  / IANDJ,KANDL,SAME
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /ROOT  / XX,U(9),W(9),NROOTS
      COMMON /SETINT/ IN(9),KN(9),NI,NJ,NK,NL,NMAX,MMAX,
     +                BP01,B00,B10,XCP00,XC00,YCP00,YC00,ZCP00,ZC00,F00,
     +                DXIJ,DYIJ,DZIJ,DXKL,DYKL,DZKL
      COMMON /SHLINF/ AG(MXGSH),CSA(MXGSH),CPA(MXGSH),CDA(MXGSH),
     *                CFA(MXGSH),CGA(MXGSH),
     *                BG(MXGSH),CSB(MXGSH),CPB(MXGSH),CDB(MXGSH),
     *                CFB(MXGSH),CGB(MXGSH),
     *                CG(MXGSH),CSC(MXGSH),CPC(MXGSH),CDC(MXGSH),
     *                CFC(MXGSH),CGC(MXGSH),
     *                DG(MXGSH),CSD(MXGSH),CPD(MXGSH),CDD(MXGSH),
     *                CFD(MXGSH),CGD(MXGSH),
     *                XI,YI,ZI,XJ,YJ,ZJ,RRI,XK,YK,ZK,XL,YL,ZL,RRK,
     *                NGA,NGB,NGC,NGD
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     +                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     +                NIJ,IJ,KL,IJKL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
C
      DIMENSION IN1(9)
C
      PARAMETER (SQRT3=1.73205080756888D+00, SQRT5=2.23606797749979D+00,
     *           SQRT7=2.64575131106459D+00, PI252=34.986836655250D+00,
     *           ZERO=0.0D+00, HALF=0.5D+00, ONE=1.0D+00)
C
C     GENERAL INTEGRAL ROUTINE FOR SPD FUNCTIONS
C
C     PI252 = PI**2.5*2
      FACTOR = PI252*QQ4
CCC      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
CCC      write(*,*)'BARIGN-start>xi,yi,zi=',xi,yi,zi
CCC      write(*,*)'BARIGN-start>xj,yj,zj=',xj,yj,zj
      NI = LIT-1
      NJ = LJT-1
      NK = 0
      NL = 0
      DXIJ = XI-XJ
      DYIJ = YI-YJ
      DZIJ = ZI-ZJ
      DXKL = ZERO
      DYKL = ZERO
      DZKL = ZERO
      NMAX = NI+NJ
      MMAX = NK+NL
      MAX = NMAX+1
      DO 100 I = 1,MAX
         N = I-1
         IF (N .LE. NI) IN1(I) = 125*N+1
         IF (N .GT. NI) IN1(I) = 125*NI+25*(N-NI)+1
  100 CONTINUE
      MAX = MMAX+1
      DO 120 K = 1,MAX
         N = K-1
         IF (N .LE. NK) KN(K) = 5*N
         IF (N .GT. NK) KN(K) = 5*NK+N-NK
  120 CONTINUE
C     KN(2) is actually accessed in XYZINT, but then never used!
      KN(2)=0
CCC      write(*,*)'BARIGN>nij,ni,nj,nmax,mmax',nij,ni,nj,nmax,mmax
CCC      write(*,*)'BARIGN>in1=',(in1(i),i=1,nmax+1)
CCC      write(*,*)'BARIGN>kn=',(kn(i),i=1,nmax+1)
C
C     ----- K,L PRIMITIVE - loops deleted
C     some variables here are maybe not needed (RRK,...)
C
      RRK = ZERO
      BRRK = ZERO
      AKXK = AK*QX
      AKYK = AK*QY
      AKZK = AK*QZ
      CSK = CKN*FACTOR
      AL = AK
      B = AK+AL
      BINV = ONE/B
      BBRRK = AL*BRRK*BINV
      CSL = CLN
      XB = QX
      YB = QY
      ZB = QZ
      BXBK = ZERO
      BYBK = ZERO
      BZBK = ZERO
      BXBI = B*(XB-XI)
      BYBI = B*(YB-YI)
      BZBI = B*(ZB-ZI)
C
C           ----- DENSITY FACTOR
C
      DKL(1) = CSK*BINV*CSL
C
C           ----- PAIR OF I,J PRIMITIVES
C
      NN = 0
      DO 440 N = 1,NIJ
         DUM = BBRRK+R(N)
CCC         IF (DUM .GT. TOL) GO TO 440
         DO 380 I = 1,IJ
            DIJ(I) = DDIJ(IJD(I)+NN)
CCC            write(*,*)'BARIGN>i,ijd(i),nn,dij(i)=',i,ijd(i),nn,dij(i)
 380     CONTINUE
         A = AA(N)
         AB = A*B
         AANDB = A+B
         EXPE = EXP(-DUM)/SQRT(AANDB)
         RHO = AB/AANDB
         XA = X1(N)
         YA = Y1(N)
         ZA = Z1(N)
         XX = RHO*((XA-XB)*(XA-XB) + (YA-YB)*(YA-YB)
     *        + (ZA-ZB)*(ZA-ZB))
CCC         write(*,*)'BARIGN>xa,ya,za,xb,yb,zb=',xa,ya,za,xb,yb,zb
CCC         write(*,*)'BARIGN>aandb,ab,a,b,dum,expe=',
CCC     $        aandb,ab,a,b,dum,expe
CCC         write(*,*)'BARIGN>rho,xx=',rho,xx
         AXAK = A*(XA-QX)
         AYAK = A*(YA-QY)
         AZAK = A*(ZA-QZ)
         AXAI = A*(XA-XI)
         AYAI = A*(YA-YI)
         AZAI = A*(ZA-ZI)
CCC         write(*,*)'BARIGN>xi,yi,zi=',xi,yi,zi
CCC         write(*,*)'BARIGN>xj,yj,zj=',xj,yj,zj
CCC         write(*,*)'BARIGN>axai,ayai,azai=',axai,ayai,azai
CCC         write(*,*)'BARIGN>bxbi,bybi,bzbi=',bxbi,bybi,bzbi
         C1X = BXBK+AXAK
         C2X = A*BXBK
         C3X = BXBI+AXAI
         C4X = B*AXAI
         C1Y = BYBK+AYAK
         C2Y = A*BYBK
         C3Y = BYBI+AYAI
         C4Y = B*AYAI
         C1Z = BZBK+AZAK
         C2Z = A*BZBK
         C3Z = BZBI+AZAI
         C4Z = B*AZAI
C
C              ----- ROOTS AND WEIGHTS FOR QUADRATURE
C
         IF (NROOTS .LE. 3) CALL RT123
         IF (NROOTS .EQ. 4) CALL ROOT4
         IF (NROOTS .EQ. 5) CALL ROOT5
         IF (NROOTS .GE. 6) CALL ROOT6
         MM = 0
         MAX = NMAX+1
C
C              COMPUTE TWO-ELECTRON INTEGRALS FOR EACH ROOT
C
         DO 420 M = 1,NROOTS
CCC            write(*,*)'BARIGN>c3x,c3y,c3z,c4x,c4y,c4z=',
CCC     $           c3x,c3y,c3z,c4x,c4y,c4z
CCC            write(*,*)'BARIGN>root,u(),w()=',m,u(m),w(m)
            U2 = U(M)*RHO
            F00 = EXPE*W(M)
            DO 400 I = 1,MAX
               IN(I) = IN1(I)+MM
 400        CONTINUE
            DUMINV = ONE/(AB+U2*AANDB)
            DM2INV = HALF*DUMINV
            BP01 = (A+U2)*DM2INV
            B00 = U2*DM2INV
            B10 = (B+U2)*DM2INV
            XCP00 = (U2*C1X+C2X)*DUMINV
            XC00 = (U2*C3X+C4X)*DUMINV
            YCP00 = (U2*C1Y+C2Y)*DUMINV
            YC00 = (U2*C3Y+C4Y)*DUMINV
            ZCP00 = (U2*C1Z+C2Z)*DUMINV
            ZC00 = (U2*C3Z+C4Z)*DUMINV
            CALL XYZINT1
            MM = MM+625
 420     CONTINUE
C
C              ----- FORM (I,J//K,L) INTEGRALS OVER FUNCTIONS
C
         CALL FORMS1(GHONDO)
 440  NN = NN+16
C
      RETURN
      END
C
      SUBROUTINE FORMS1(GHONDO)
C-----------------------------------------------------------------------
C
C     This routine is taken from int2a.src(forms) and modified for blur
C     functions. Note that IJX(),IJY(),IJZ() are different than
C     the ones in SUBROUTINE BLURIN, because XIN(),YIN(),ZIN()
C     are calculated that way....
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C
      DIMENSION GHONDO(*)
C
      COMMON /DENS  / DKL(225),DIJ(225)
      COMMON /INTDEX/ IJGT(225),IJX(225),IJY(225),IJZ(225),IK(225),
     *                KLGT(225),KLX(225),KLY(225),KLZ(225)
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJ,IJ,KL,IJKL
      COMMON /ROOT  / XX,U(9),W(9),NROOTS
      COMMON /XYZ   / XIN(5625),YIN(5625),ZIN(5625)
C
C     ----- FORM INTEGRALS OVER FUNCTIONS -----
C     DIMENSIONING XIN(81,5), AND ROLLING UP THE COMPUTATION
C     OF GHONDO IN A LOOP OF LENGTH NROOTS ADDS 33 SECONDS TO
C     A 240 SECOND INTEGRAL COMPUTATION JOB.  LEAVE IT UNROLLED.
C
CCC      write(*,*)'FORMS1>nroots=',nroots
      GO TO (100,140,180,220,260,300,340,380,420),NROOTS
C
  100 DO 120 I = 1,IJ
      D1 = DIJ(I)
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
CCC      write(*,*)'FORMS1>i,nx,ny,nz,d1,dkl,ghondo(i),xyzINT=',
CCC     $     i,nx,ny,nz,d1,dkl(1),ghondo(i),
CCC     $     ( XIN(NX )*YIN(NY )*ZIN(NZ ) )*D1*DKL(1)
  120 GHONDO(I)=( XIN(NX )*YIN(NY )*ZIN(NZ ) )*D1*DKL(1)+GHONDO(I)
      RETURN
C
  140 DO 160 I = 1,IJ
      D1 = DIJ(I)
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
CCC      write(*,*)'FORMS1>i,nx,ny,nz,d1,dkl,ghondo(i),xyzINT=',
CCC     $     i,nx,ny,nz,d1,dkl(1),ghondo(i),
CCC     $     ( XIN(NX )*YIN(NY )*ZIN(NZ ) +XIN(NX+625)*YIN(NY+625)*
CCC     +     ZIN(NZ+625) )*D1*DKL(1)
  160 GHONDO(I)=( XIN(NX )*YIN(NY )*ZIN(NZ ) +XIN(NX+625)*YIN(NY+625)*
     +     ZIN(NZ+625) )*D1*DKL(1)+GHONDO(I)
      RETURN
C
  180 DO 200 I = 1,IJ
      D1 = DIJ(I)
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
  200 GHONDO(I) = ( XIN(NX )*YIN(NY )*ZIN(NZ ) +XIN(NX+625)*YIN(NY+625)*
     +     ZIN(NZ+625) +XIN(NX+1250)*YIN(NY+1250)*ZIN(NZ+1250) )*D1*
     +     DKL(1)+GHONDO(I)
      RETURN
C
  220 DO 240 I = 1,IJ
      D1 = DIJ(I)
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
  240 GHONDO(I) = ( XIN(NX )*YIN(NY )*ZIN(NZ ) +XIN(NX+625)*YIN(NY+625)*
     +     ZIN(NZ+625) +XIN(NX+1250)*YIN(NY+1250)*ZIN(NZ+1250) +XIN(NX+
     +     1875)*YIN(NY+1875)*ZIN(NZ+1875) )*D1*DKL(1)+GHONDO(I)
      RETURN
C
  260 DO 280 I = 1,IJ
      D1 = DIJ(I)
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
  280 GHONDO(I) = ( XIN(NX )*YIN(NY )*ZIN(NZ ) +XIN(NX+625)*YIN(NY+625)*
     +     ZIN(NZ+625) +XIN(NX+1250)*YIN(NY+1250)*ZIN(NZ+1250) +XIN(NX+
     +     1875)*YIN(NY+1875)*ZIN(NZ+1875) +XIN(NX+2500)*YIN(NY+2500)*
     +     ZIN(NZ+2500) )*D1*DKL(1)+GHONDO(I)
      RETURN
C
  300 DO 320 I = 1,IJ
      D1 = DIJ(I)
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
  320 GHONDO(I) = ( XIN(NX )*YIN(NY )*ZIN(NZ ) +XIN(NX+625)*YIN(NY+625)*
     +     ZIN(NZ+625) +XIN(NX+1250)*YIN(NY+1250)*ZIN(NZ+1250) +XIN(NX+
     +     1875)*YIN(NY+1875)*ZIN(NZ+1875) +XIN(NX+2500)*YIN(NY+2500)*
     +     ZIN(NZ+2500) +XIN(NX+3125)*YIN(NY+3125)*ZIN(NZ+3125))
     +     *D1*DKL(1)+GHONDO(I)
      RETURN
C
  340 DO 360 I = 1,IJ
      D1 = DIJ(I)
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
  360 GHONDO(I) = ( XIN(NX )*YIN(NY )*ZIN(NZ ) +XIN(NX+625)*YIN(NY+625)*
     +     ZIN(NZ+625) +XIN(NX+1250)*YIN(NY+1250)*ZIN(NZ+1250) +XIN(NX+
     +     1875)*YIN(NY+1875)*ZIN(NZ+1875) +XIN(NX+2500)*YIN(NY+2500)*
     +     ZIN(NZ+2500) +XIN(NX+3125)*YIN(NY+3125)*ZIN(NZ+3125) +XIN(NX+
     +     3750)*YIN(NY+3750)*ZIN(NZ+3750))
     +     *D1*DKL(1)+GHONDO(I)
      RETURN
C
  380 DO 400 I = 1,IJ
      D1 = DIJ(I)
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
  400 GHONDO(I) = ( XIN(NX )*YIN(NY )*ZIN(NZ ) +XIN(NX+625)*YIN(NY+625)*
     +     ZIN(NZ+625) +XIN(NX+1250)*YIN(NY+1250)*ZIN(NZ+1250) +XIN(NX+
     +     1875)*YIN(NY+1875)*ZIN(NZ+1875) +XIN(NX+2500)*YIN(NY+2500)*
     +     ZIN(NZ+2500) +XIN(NX+3125)*YIN(NY+3125)*ZIN(NZ+3125) +XIN(NX+
     +     3750)*YIN(NY+3750)*ZIN(NZ+3750) +XIN(NX+4375)*YIN(NY+4375)*
     +     ZIN(NZ+4375))*D1*DKL(1)+GHONDO(I)
      RETURN
C
  420 DO 440 I = 1,IJ
      D1 = DIJ(I)
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
  440 GHONDO(I) = ( XIN(NX )*YIN(NY )*ZIN(NZ ) +XIN(NX+625)*YIN(NY+625)*
     +     ZIN(NZ+625) +XIN(NX+1250)*YIN(NY+1250)*ZIN(NZ+1250) +XIN(NX+
     +     1875)*YIN(NY+1875)*ZIN(NZ+1875) +XIN(NX+2500)*YIN(NY+2500)*
     +     ZIN(NZ+2500) +XIN(NX+3125)*YIN(NY+3125)*ZIN(NZ+3125) +XIN(NX+
     +     3750)*YIN(NY+3750)*ZIN(NZ+3750) +XIN(NX+4375)*YIN(NY+4375)*
     +     ZIN(NZ+4375) +XIN(NX+5000)*YIN(NY+5000)*ZIN(NZ+5000))
     +     *D1*DKL(1)+GHONDO(I)
      RETURN
      END
C
      SUBROUTINE IJPRIM1(DDIJ)
C------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C
      LOGICAL IANDJ,KANDL,SAME,OUT,NORM
C
      PARAMETER (MXGSH=30, MXG2=MXGSH*MXGSH)
C
      DIMENSION DDIJ(16*MXG2)
C
      COMMON /IJGNRL/ A(MXG2),R(MXG2),X1(MXG2),Y1(MXG2),Z1(MXG2),
     *                IJD(225)
      COMMON /MISC  / IANDJ,KANDL,SAME
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /SHLINF/ AG(MXGSH),CSA(MXGSH),CPA(MXGSH),CDA(MXGSH),
     *                CFA(MXGSH),CGA(MXGSH),
     *                BG(MXGSH),CSB(MXGSH),CPB(MXGSH),CDB(MXGSH),
     *                CFB(MXGSH),CGB(MXGSH),
     *                CG(MXGSH),CSC(MXGSH),CPC(MXGSH),CDC(MXGSH),
     *                CFC(MXGSH),CGC(MXGSH),
     *                DG(MXGSH),CSD(MXGSH),CPD(MXGSH),CDD(MXGSH),
     *                CFD(MXGSH),CGD(MXGSH),
     *                XI,YI,ZI,XJ,YJ,ZJ,RRI,XK,YK,ZK,XL,YL,ZL,RRK,
     *                NGA,NGB,NGC,NGD
       COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJ,IJ,KL,IJKL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
C
      PARAMETER (SQRT3=1.73205080756888D+00,
     *           SQRT5=2.23606797749979D+00,
     *           SQRT7=2.64575131106459D+00,
     *           ZERO=0.0D+00, ONE=1.0D+00)
C
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
      TOL=75.0D+00
      MAX = MAXJ
      N = 0
      NN = 0
      NM = -2**20
      DO 180 I = MINI,MAXI
         GO TO (100,100,120,120,100,120,120,100,120,120,
     1          100,120,120,100,120,120,120,120,120,100,
     1          100,120,120,100,120,120,120,120,120,100,
     1          120,120,100,120,120),I
  100    NM = NN
  120    NN = NM
         IF (IANDJ) MAX = I
         DO 170 J = MINJ,MAX
            GO TO (140,140,160,160,140,160,160,140,160,160,
     1             140,160,160,140,160,160,160,160,160,140,
     1             140,160,160,140,160,160,160,160,160,140,
     1             160,160,140,160,160),J
  140       NN = NN+1
  160       N = N+1
            IJD(N) = NN
  170    CONTINUE
  180 CONTINUE
CCC      write(*,*)'IJPRIM1>nij,nn=',nij,nn
CCC      write(*,*)'IJPRIM1>mini,maxi,minj,maxj=',
CCC     $     mini,maxi,minj,maxj
CCC      write(*,*)'IJPRIM1>mink,maxk,minl,maxl=',
CCC     $     mink,maxk,minl,maxl
CCC      write(*,*)'IJPRIM1>nga,ngb=',nga,ngb
CCC      write(*,*)'IJPRIM1>n,ijd=',n,(ijd(i),i=1,n)
C
C     ----- I PRIMITIVE
C
      NIJ = 0
      JBMAX = NGB
      DO 540 IA = 1,NGA
         AI = AG(IA)
         ARRI = AI*RRI
         AXI = AI*XI
         AYI = AI*YI
         AZI = AI*ZI
         CSI = CSA(IA)
         CPI = CPA(IA)
         CDI = CDA(IA)
         CFI = CFA(IA)
         CGI = CGA(IA)
C
C        ----- J PRIMITIVE
C
         IF (IANDJ) JBMAX = IA
         DO 520 JB = 1,JBMAX
            AJ = BG(JB)
            AA = AI+AJ
            AAINV = ONE/AA
            DUM = AJ*ARRI*AAINV
            IF (DUM .GT. TOL) GO TO 520
            CSJ = CSB(JB)
            CPJ = CPB(JB)
            CDJ = CDB(JB)
            CFJ = CFB(JB)
            CGJ = CGB(JB)
            NM = 16*NIJ
            NN = NM
            NIJ = NIJ+1
            R(NIJ) = DUM
            A(NIJ) = AA
            X1(NIJ) = (AXI+AJ*XJ)*AAINV
            Y1(NIJ) = (AYI+AJ*YJ)*AAINV
            Z1(NIJ) = (AZI+AJ*ZJ)*AAINV
C
C           ----- DENSITY FACTOR
C
            DUM1 = ZERO
            DUM2 = ZERO
            DO 420 I = MINI,MAXI
               GO TO (200,220,420,420,240,420,420,260,420,420,
     1                261,420,420,262,420,420,420,420,420,263,
     1                264,420,420,265,420,420,420,420,420,266,
     1                420,420,267,420,420),I
  200          DUM1 = CSI*AAINV
               GO TO 280
  220          DUM1 = CPI*AAINV
               GO TO 280
  240          DUM1 = CDI*AAINV
               GO TO 280
  260          IF (NORM) DUM1 = DUM1*SQRT3
               GO TO 280
  261          DUM1 = CFI*AAINV
               GO TO 280
  262          IF (NORM) DUM1 = DUM1*SQRT5
               GO TO 280
  263          IF (NORM) DUM1 = DUM1*SQRT3
               GO TO 280
  264          DUM1 = CGI*AAINV
               GO TO 280
  265          IF (NORM) DUM1 = DUM1*SQRT7
               GO TO 280
  266          IF (NORM) DUM1 = DUM1*SQRT5/SQRT3
               GO TO 280
  267          IF (NORM) DUM1 = DUM1*SQRT3
  280          IF (IANDJ) MAX = I
               DO 400 J = MINJ,MAX
                  GO TO (300,320,400,400,340,400,400,360,400,400,
     1                   361,400,400,362,400,400,400,400,400,363,
     1                   364,400,400,365,400,400,400,400,400,366,
     1                   400,400,367,400,400),J
  300             DUM2 = DUM1*CSJ
                  GO TO 380
  320             DUM2 = DUM1*CPJ
                  GO TO 380
  340             DUM2 = DUM1*CDJ
                  GO TO 380
  360             IF (NORM) DUM2 = DUM2*SQRT3
                  GO TO 380
  361             DUM2 = DUM1*CFJ
                  GO TO 380
  362             IF (NORM) DUM2 = DUM2*SQRT5
                  GO TO 380
  363             IF (NORM) DUM2 = DUM2*SQRT3
                  GO TO 380
  364             DUM2 = DUM1*CGJ
                  GO TO 380
  365             IF (NORM) DUM2 = DUM2*SQRT7
                  GO TO 380
  366             IF (NORM) DUM2 = DUM2*SQRT5/SQRT3
                  GO TO 380
  367             IF (NORM) DUM2 = DUM2*SQRT3
  380             NN = NN+1
                  DDIJ(NN) = DUM2
CCC                  write(*,*)'IJPRIM1>aa,nn,ddij(nn)=',aa,nn,ddij(nn)
CCC                  write(*,*)'IJPRIM1>csi,csj,cpi,cpj=',csi,csj,cpi,cpj
  400          CONTINUE
  420       CONTINUE
CCC            write(*,*)'IPRIM1>iandj,ai,jb,lit,mini',
CCC     $           iandj,ai,jb,lit,mini
            IF ( .NOT. IANDJ) GO TO 520
            IF (IA .EQ. JB) GO TO 520
            GO TO (500,440,460,455,450),LIT
  440       IF (MINI .EQ. 2) GO TO 500
            DDIJ(NM+2) = DDIJ(NM+2)+CSI*CPJ*AAINV
CCC            write(*,*)'IJPRIM1>nm+2,ddij(nm+2)=',nm+2,ddij(nm+2)
            GO TO 480
  450       DDIJ(NM+10) = DDIJ(NM+10)+DDIJ(NM+10)
            DDIJ(NM+9) = DDIJ(NM+9)+DDIJ(NM+9)
            DDIJ(NM+8) = DDIJ(NM+8)+DDIJ(NM+8)
            DDIJ(NM+7) = DDIJ(NM+7)+DDIJ(NM+7)
  455       DDIJ(NM+6) = DDIJ(NM+6)+DDIJ(NM+6)
            DDIJ(NM+5) = DDIJ(NM+5)+DDIJ(NM+5)
            DDIJ(NM+4) = DDIJ(NM+4)+DDIJ(NM+4)
  460       DDIJ(NM+2) = DDIJ(NM+2)+DDIJ(NM+2)
  480       DDIJ(NM+3) = DDIJ(NM+3)+DDIJ(NM+3)
  500       DDIJ(NM+1) = DDIJ(NM+1)+DDIJ(NM+1)
CCC            write(*,*)'IJPRIM1>nm+10,ddij(nm+10)=',nm+2,ddij(nm+10)
CCC            write(*,*)'IJPRIM1>nm+9,ddij(nm+9)=',nm+9,ddij(nm+9)
CCC            write(*,*)'IJPRIM1>nm+8,ddij(nm+8)=',nm+8,ddij(nm+8)
CCC            write(*,*)'IJPRIM1>nm+7,ddij(nm+7)=',nm+7,ddij(nm+7)
CCC            write(*,*)'IJPRIM1>nm+6,ddij(nm+6)=',nm+6,ddij(nm+6)
CCC            write(*,*)'IJPRIM1>nm+5,ddij(nm+5)=',nm+5,ddij(nm+5)
CCC            write(*,*)'IJPRIM1>nm+4,ddij(nm+4)=',nm+4,ddij(nm+4)
CCC            write(*,*)'IJPRIM1>nm+2,ddij(nm+2)=',nm+2,ddij(nm+2)
CCC            write(*,*)'IJPRIM1>nm+3,ddij(nm+3)=',nm+3,ddij(nm+3)
CCC            write(*,*)'IJPRIM1>nm+1,ddij(nm+1)=',nm+1,ddij(nm+1)
  520    CONTINUE
  540 CONTINUE
CCC      write(*,*)'IJPRIM1>(end):nij=',nij
      RETURN
      END
C
      SUBROUTINE XYZINT1
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C
      LOGICAL N0,N1,M0,M1
C
      COMMON /SETINT/ I(9),K(9),NIMAX,NJMAX,NKMAX,NLMAX,NMAX,MMAX
     +               ,BP01,B00,B10,XCP00,XC00,YCP00,YC00,ZCP00,ZC00,F00
     +               ,DXIJ,DYIJ,DZIJ,DXKL,DYKL,DZKL
      COMMON /XYZ   / XINT(5625),YINT(5625),ZINT(5625)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      N0 = NMAX .EQ. 0
      N1 = NMAX .LE. 1
      M0 = MMAX .EQ. 0
      M1 = MMAX .LE. 1
CCC      write(*,*)'XYZINT1>n0,n1,m0,m1=',n0,n1,m0,m1
CCC      write(*,*)'XYZINT1>i=',(i(j),j=1,9)
CCC      write(*,*)'XYZINT1>k=',(k(j),j=1,9)
CCC      write(*,*)'XYZINT1>NIMAX,NJMAX,NKMAX,NLMAX,NMAX,MMAX=',
CCC     $     NIMAX,NJMAX,NKMAX,NLMAX,NMAX,MMAX
CCC      write(*,*)'XYZINT1>BP01,B00,B10,XCP00,XC00,YCP00,YC00,ZCP00,ZC00,
CCC     $     F00=',BP01,B00,B10,XCP00,XC00,YCP00,YC00,ZCP00,ZC00,F00
CCC       write(*,*)'XYZINT1>DXIJ,DYIJ,DZIJ,DXKL,DYKL,DZKL=',
CCC     $     DXIJ,DYIJ,DZIJ,DXKL,DYKL,DZKL
C
C     ----- I(0,0) -----
C
      I1 = I(1)
      XINT(I1) = ONE
      YINT(I1) = ONE
      ZINT(I1) = F00
CCC      write(*,*)'XYZINT1>i1,x,y,z=',i1,xint(i1),yint(i1),zint(i1)
      IF (N0 .AND. M0) RETURN
      I2 = I(2)
      K2 = K(2)
      CP10 = B00
      IF (N0) GO TO 100
C
C     ----- I(1,0) -----
C
      XINT(I2) = XC00
      YINT(I2) = YC00
      ZINT(I2) = ZC00*F00
CCC      write(*,*)'XYZINT1>i2,x,y,z=',i2,xint(i2),yint(i2),zint(i2)
      IF (M0) GO TO 120
C
C     ----- I(0,1) -----
C
  100 I3 = I1+K2
      XINT(I3) = XCP00
      YINT(I3) = YCP00
      ZINT(I3) = ZCP00*F00
      IF (N0) GO TO 120
C
C     ----- I(1,1) -----
C
      I3 = I2+K2
      XINT(I3) = XCP00*XINT(I2)+CP10
      YINT(I3) = YCP00*YINT(I2)+CP10
      ZINT(I3) = ZCP00*ZINT(I2)+CP10*F00
  120 IF (N1) GO TO 180
      C10 = ZERO
      I3 = I1
      I4 = I2
      DO 160 N = 2,NMAX
      C10 = C10+B10
C
C     ----- I(N,0) -----
C
      I5 = I(N+1)
      XINT(I5) = C10*XINT(I3)+XC00*XINT(I4)
      YINT(I5) = C10*YINT(I3)+YC00*YINT(I4)
      ZINT(I5) = C10*ZINT(I3)+ZC00*ZINT(I4)
CCC      write(*,*)'XYZINT1>i5,x,y,z=',i5,xint(i5),yint(i5),zint(i5)
      IF (M0) GO TO 140
      CP10 = CP10+B00
C
C     ----- I(N,1) -----
C
      I3 = I5+K2
      XINT(I3) = XCP00*XINT(I5)+CP10*XINT(I4)
      YINT(I3) = YCP00*YINT(I5)+CP10*YINT(I4)
      ZINT(I3) = ZCP00*ZINT(I5)+CP10*ZINT(I4)
  140 I3 = I4
  160 I4 = I5
  180 IF (M1) GO TO 240
      CP01 = ZERO
      C01 = B00
      I3 = I1
      I4 = I1+K2
      DO 220 M = 2,MMAX
      CP01 = CP01+BP01
C
C     ----- I(0,M) -----
C
      I5 = I1+K(M+1)
      XINT(I5) = CP01*XINT(I3)+XCP00*XINT(I4)
      YINT(I5) = CP01*YINT(I3)+YCP00*YINT(I4)
      ZINT(I5) = CP01*ZINT(I3)+ZCP00*ZINT(I4)
      IF (N0) GO TO 200
      C01 = C01+B00
C
C     ----- I(1,M) -----
C
      I3 = I2+K(M+1)
      XINT(I3) = XC00*XINT(I5)+C01*XINT(I4)
      YINT(I3) = YC00*YINT(I5)+C01*YINT(I4)
      ZINT(I3) = ZC00*ZINT(I5)+C01*ZINT(I4)
  200 I3 = I4
  220 I4 = I5
  240 IF (N1 .OR. M1) GO TO 300
C
C     ----- I(N,M) -----
C
      C01 = B00
      K3 = K2
      DO 280 M = 2,MMAX
      K4 = K(M+1)
      C01 = C01+B00
      I3 = I1
      I4 = I2
      C10 = B10
      DO 260 N = 2,NMAX
      I5 = I(N+1)
      XINT(I5+K4) = C10*XINT(I3+K4)+XC00*XINT(I4+K4)+C01*XINT(I4+K3)
      YINT(I5+K4) = C10*YINT(I3+K4)+YC00*YINT(I4+K4)+C01*YINT(I4+K3)
      ZINT(I5+K4) = C10*ZINT(I3+K4)+ZC00*ZINT(I4+K4)+C01*ZINT(I4+K3)
      C10 = C10+B10
      I3 = I4
  260 I4 = I5
  280 K3 = K4
  300 IF (NJMAX .EQ. 0) GO TO 440
C
C     ----- I(NI,NJ,M) -----
C
      M = 0
      I5 = I(NMAX+1)
  320 MIN = NIMAX
      KM = K(M+1)
  340 N = NMAX
      I3 = I5+KM
  360 I4 = I(N)+KM
      XINT(I3) = XINT(I3)+DXIJ*XINT(I4)
      YINT(I3) = YINT(I3)+DYIJ*YINT(I4)
      ZINT(I3) = ZINT(I3)+DZIJ*ZINT(I4)
CCC      write(*,*)'XYZINT1>i3,x,y,z=',i3,xint(i3),yint(i3),zint(i3)
      I3 = I4
      N = N-1
      IF (N .GT. MIN) GO TO 360
      MIN = MIN+1
      IF (MIN .LT. NMAX) GO TO 340
      IF (NIMAX .EQ. 0) GO TO 420
      I3 = 25+KM+I1
      DO 400 NJ = 1,NJMAX
      I4 = I3
      DO 380 NI = 1,NIMAX
      XINT(I4) = XINT(I4+100)+DXIJ*XINT(I4-25)
      YINT(I4) = YINT(I4+100)+DYIJ*YINT(I4-25)
      ZINT(I4) = ZINT(I4+100)+DZIJ*ZINT(I4-25)
CCC      write(*,*)'XYZINT1>i4,x,y,z=',i4,xint(i4),yint(i4),zint(i4)
  380 I4 = I4+125
  400 I3 = I3+25
  420 M = M+1
      IF (M .LE. MMAX) GO TO 320
  440 IF (NLMAX .EQ. 0) GO TO 600
C
C     ----- I(NI,NJ,NK,NL) -----
C
      I5 = K(MMAX+1)
      IA = I1
      NI = 0
  460 NJ = 0
      IB = IA
  480 MIN = NKMAX
  500 M = MMAX
      I3 = IB+I5
  520 I4 = IB+K(M)
      XINT(I3) = XINT(I3)+DXKL*XINT(I4)
      YINT(I3) = YINT(I3)+DYKL*YINT(I4)
      ZINT(I3) = ZINT(I3)+DZKL*ZINT(I4)
      I3 = I4
      M = M-1
      IF (M .GT. MIN) GO TO 520
      MIN = MIN+1
      IF (MIN .LT. MMAX) GO TO 500
      IF (NKMAX .EQ. 0) GO TO 580
      I3 = IB+1
      DO 560 NL = 1,NLMAX
      I4 = I3
      DO 540 NK = 1,NKMAX
      XINT(I4) = XINT(I4+4)+DXKL*XINT(I4-1)
      YINT(I4) = YINT(I4+4)+DYKL*YINT(I4-1)
      ZINT(I4) = ZINT(I4+4)+DZKL*ZINT(I4-1)
  540 I4 = I4+5
  560 I3 = I3+1
  580 NJ = NJ+1
      IB = IB+25
      IF (NJ .LE. NJMAX) GO TO 480
      NI = NI+1
      IA = IA+125
      IF (NI .LE. NIMAX) GO TO 460
  600 CONTINUE
      RETURN
C
##ELSE
      SUBROUTINE BLUR_NULL
      RETURN
##ENDIF
C
      END
