C 11 DEC 06 - MWS - UPDATE COMMON BLOCKS FOR FMO CHANGES
C 19 NOV 06 - PEA - NJANDK: GENERATE 4 CENTER OVERLAP INTS IF ASKED
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 25 OCT 06 - SPW - CHANGES FOR NEO VERSION 2
C 25 MAR 06 - SPW - NEW MODULE FOR NEO INTEGRALS
C
C*MODULE NEOINT  *DECK NUCHST
      SUBROUTINE NUCHST(H,S,T,Z,ESP1E,LL2,SOME,LZINT,DOESP,DBUG)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION NINE,MOROKM
C
      LOGICAL SOME,DBUG,LZINT,DOESP
      LOGICAL IANDJ,DOUBLE,GOPARR,DSKWRK,MASWRK,SCREEN
C
C     COSMO CHANGES
C
      LOGICAL ISEPS, USEPS
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MAXSH=84)
C
      DIMENSION H(LL2),S(LL2),T(LL2),Z(LL2),ESP1E(LL2)
      DIMENSION SBLK(784),TBLK(784),VBLK(784),ZBLK(784),
     *          FT(784),DIJ(784),IJX(784),IJY(784),IJZ(784),
     *          XIN(343),YIN(343),ZIN(343),CONI(MAXSH),CONJ(MAXSH),
     *          IX(84),IY(84),IZ(84),JX(84),JY(84),JZ(84)
C
      PARAMETER (MXCHRM=1)
C
      COMMON /CHMGMS/ XCHM(MXCHRM),YCHM(MXCHRM),ZCHM(MXCHRM),
     *                DXELMM(MXCHRM),DYELMM(MXCHRM),DZELMM(MXCHRM),
     *                QCHM(MXCHRM),NCHMAT,KCHRMM
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NEOSHL/ NGAUSS,NGAUSE,NGAUSN,NTSHEL,NNUCSH
      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 /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      COMMON /NUCMOR/ QMNUCM(20)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCINP / VLAMB,SCREEN
      COMMON /STV   / XINT,YINT,ZINT,TAA,X0,Y0,Z0,
     *                XI,YI,ZI,XJ,YJ,ZJ,NI,NJ
      COMMON /SHLNRM/ PNRM(84)
      COMMON /SYMIND/ TOL,II,JJ,LIT,LJT,MINI,MINJ,MAXI,MAXJ,IANDJ
      COMMON /WERMIT/ W1,W2(2),W3(3),W4(4),W5(5),W6(6),W7(7),
     *                W8(8),W9(9),W10(10)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
C     COSMO CHANGES
C
      PARAMETER (MAXDEN=25*MXATM, LENABC=2000, NPPA=1082)
      COMMON /CONV  / DENTOL,EN,ETOT,EHF,EHF0,DIFF,ITER,ICALP,ICBET
      COMMON /COSMO3/ COSZAN(NPPA),CORZAN(3,NPPA)
      COMMON /ISEPS / ISEPS,USEPS
      COMMON /SOLV  / FEPSI,RDS,DISEX2,COSURF(3,LENABC),SRAD(MXATM),
     *                QDEN(MAXDEN),AR(LENABC),
     *                NSPA,NPS,NPS2,NDEN,NPSPHER
c
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /FMOPNT/ LICHFG,LMULFG,LIDMREC,LFRGNAM,LLAYFRG,LINDAT,
     *                LNCBS,LFMOZAN,LFMOC,LFMOMAS,LIZBAS,LIAGLOB,LIBDGH,
     *                LIABDFG,LJABDFG,LNCAO,LIDXCAO,LIAPRJO,LJAPRJO,
     *                LCOREAO,LOCCCOR,LSHIFTB,LIODFMO,LFMODA,LFMODB,
     *                LFMOESPA,LFMOESPB,LLOCFMO,LSCFFRG,LRIJ,LPOPMUL,
     *                LPOPMAT,LIALOC,LINDBD,LIATFRG,LINDFRG,LINDGFRG,
     *                LNATFRG,LNAT0FRG,LIANFRG,LZANFRG,LCFRG,LLIBISH,
     *                LLIBNSH,LLIBNG,LINDATG,LFMOBUF(3),LFMODE,LNUMFRG,
     *                LLOCTAT,LIAOGLOB,LLOADM,LFMOGE,LDGRID,LIODCFMO,
     *                LJOB2GRP,LFMOPG,LEMOCDR,LUNTXYZ,LUNTROT,LSTONEP,
     *                LMAPSU,LFRGMUL
      COMMON /FMORUN/ ESPSCF,E0SCF(2),IDAFMO,ICURLAY,ICURFG,JCURFG,
     *                KCURFG,NCURSH,NGAU,ICURPOP,IFMOSTP,MONCOR,NAT1E,
     *                NEEDR,MODRST,NORBPROJ,NUNESP,ISKIPESP,ICURUNT
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,
     *           THIRTN=13.0D+00, FIFTEN=15.0D+00,
     *           PI212=1.1283791670955D+00, RLN10=2.30258D+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,
     *          5, 0, 0, 4, 4, 1, 0, 1, 0, 3,
     *          3, 2, 0, 2, 0, 3, 1, 1, 2, 2,
     *          1,
     *          6, 0, 0, 5, 5, 1, 0, 1, 0, 4,
     *          4, 2, 0, 2, 0, 4, 1, 1, 3, 3,
     *          0, 3, 3, 2, 1, 2, 1, 2/
      DATA IX / 1, 8, 1, 1,15, 1, 1, 8, 8, 1,
     *         22, 1, 1,15,15, 8, 1, 8, 1, 8,
     *         29, 1, 1,22,22, 8, 1, 8, 1,15,
     *         15, 1,15, 8, 8,
     *         36, 1, 1,29,29, 8, 1, 8, 1,22,
     *         22,15, 1,15, 1,22, 8, 8,15,15,
     *          8,
     *         43, 1, 1,36,36, 8, 1, 8, 1,29,
     *         29,15, 1,15, 1,29, 8, 8,22,22,
     *          1,22,22,15, 8,15, 8,15/
      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,
     *          0, 5, 0, 1, 0, 4, 4, 0, 1, 2,
     *          0, 3, 3, 0, 2, 1, 3, 1, 2, 1,
     *          2,
     *          0, 6, 0, 1, 0, 5, 5, 0, 1, 2,
     *          0, 4, 4, 0, 2, 1, 4, 1, 3, 0,
     *          3, 2, 1, 3, 3, 1, 2, 2/
      DATA IY / 1, 1, 8, 1, 1,15, 1, 8, 1, 8,
     *          1,22, 1, 8, 1,15,15, 1, 8, 8,
     *          1,29, 1, 8, 1,22,22, 1, 8,15,
     *          1,15, 8,15, 8,
     *          1,36, 1, 8, 1,29,29, 1, 8,15,
     *          1,22,22, 1,15, 8,22, 8,15, 8,
     *         15,
     *          1,43, 1, 8, 1,36,36, 1, 8,15,
     *          1,29,29, 1,15, 8,29, 8,22, 1,
     *         22,15, 8,22,22, 8,15,15/
      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,
     *          0, 0, 5, 0, 1, 0, 1, 4, 4, 0,
     *          2, 0, 2, 3, 3, 1, 1, 3, 1, 2,
     *          2,
     *          0, 0, 6, 0, 1, 0, 1, 5, 5, 0,
     *          2, 0, 2, 4, 4, 1, 1, 4, 0, 3,
     *          3, 1, 2, 1, 2, 3, 3, 2/
      DATA IZ / 1, 1, 1, 8, 1, 1,15, 1, 8, 8,
     *          1, 1,22, 1, 8, 1, 8,15,15, 8,
     *          1, 1,29, 1, 8, 1, 8,22,22, 1,
     *         15,15, 8, 8,15,
     *          1, 1,36, 1, 8, 1, 8,29,29, 1,
     *         15, 1,15,22,22, 8, 8,22, 8,15,
     *         15,
     *          1, 1,43, 1, 8, 1, 8,36,36, 1,
     *         15, 1,15,29,29, 8, 8,29, 1,22,
     *         22, 8,15, 8,15,22,22,15/
C
      DATA MOROKM/8HMOROKUMA/
      DATA NONE/4HNONE/
C
C     ----- COMPUTE CONVENTIONAL H, S, AND T INTEGRALS -----
C
      TOL = RLN10*ITOL
C     NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
      IF(SOME  .AND.  SCREEN) WRITE(IW,9010) VLAMB
      IF(DOESP.AND.LZINT) THEN
C        ZBLK IS REUSED FOR SOME OTHER PURPOSE.
         WRITE(IW,*) 'STORAGE NOT ALLOCATED IN HSANDT'
         CALL ABRT
      ENDIF
C
C     ----- MOPAC INTEGRALS ARE DONE ELSEWHERE -----
C
      IF(MPCTYP.NE.NONE) THEN
         CALL MPCINT
         RETURN
      END IF
C     COSMO CHANGES VOLKER FOR DIRECT SCF, MRZ 99
C
C     FOR DIRECT SCF, WE NEED THE DIFFERENCE OF H BETWEEN THE SCF
C     ITERATIONS IN RHFCL, SO HOLD IS COPIED INTO SECTION 87
C     AND HNEW IS IN SECTION 11 AS USUAL, BEGINNING IN SCF-CYCLE 2
C
      IF(ISEPS  .AND.  ITER.GT.0) THEN
         CALL DAREAD(IDAF,IODA,H,LL2,11,0)
         CALL DAWRIT(IDAF,IODA,H,LL2,87,0)
      END IF
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 = NSHELL + 1
         IEND   = NTSHEL
         JSTART = NSHELL + 1
         LOCIJ  = 0
         NATST  = 1
         NATED  = NAT+NCHMAT
         IF(ISEPS) NATED = NAT+NPS
         IF(DOESP) NATED=NAT1E
         ISAVE  = 0
         L1 = NUMNB
         L2 = (NUMNB*(NUMNB+1))/2
      END IF
C
C     ---- NEO ---
C     ---- SCALE KINETIC ENERGY BY MASS OF NUCLEUS RELATIVE TO E- ----
C
      SCAL = ONE/QMNUCM(1)
C
      CALL VCLR(H,1,LL2)
      CALL VCLR(S,1,LL2)
      CALL VCLR(T,1,LL2)
      IF(LZINT) CALL VCLR(Z,1,LL2)
      IF(DOESP) CALL VCLR(ESP1E,1,LL2)
C
C     ----- INTIALIZE PARALLEL -----
C
      IPCOUNT = ME - 1
C
C     ----- I SHELL -----
C
      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
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
            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.GE. 2).AND.(J.LE. 4)) FT(IJ) = FIVE
                  IF ((J.GE. 5).AND.(J.LE.10)) FT(IJ) = SEVEN
                  IF ((J.GE.11).AND.(J.LE.20)) FT(IJ) = NINE
                  IF ((J.GE.21).AND.(J.LE.35)) FT(IJ) = ELEVEN
                  IF ((J.GE.36).AND.(J.LE.56)) FT(IJ) = THIRTN
                  IF ((J.GE.57).AND.(J.LE.84)) FT(IJ) = FIFTEN
  140          CONTINUE
  160       CONTINUE
C
            CALL VCLR(SBLK,1,IJ)
            CALL VCLR(TBLK,1,IJ)
            CALL VCLR(VBLK,1,IJ)
            IF(LZINT.OR.DOESP) CALL VCLR(ZBLK,1,IJ)
C
C     ----- I PRIMITIVE
C
            JGMAX = J2
            DO 520 IG = I1,I2
               AI = EX(IG)
               ARRI = AI*RR
               AXI = AI*XI
               AYI = AI*YI
               AZI = AI*ZI
               CALL SETCONI(CONI,IG)
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
                  IF (DUM .GT. TOL) GO TO 500
                  FAC = EXP(-DUM)
                  CALL SETCONI(CONJ,JG)
                  AX = (AXI+AJ*XJ)*AA1
                  AY = (AYI+AJ*YJ)*AA1
                  AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
C     THE ONLY REASON WHY -ILZ WORKS WITH THIS DENSITY THAT ASSUMES
C     HERMITICITY IS BECAUSE <I|-ILZ|I>=0 (MOMENTUM QUENCHING).
C
                  DOUBLE=IANDJ.AND.IG.NE.JG
                  MAX = MAXJ
                  NN = 0
                  DTWO=ONE
                  IF(DOUBLE) DTWO=TWO
C                 NOTE THAT PNRM FACTORS FOR S AND P SHELLS ARE 1.
                  SPDIJ=CS(IG)*CP(JG)*FAC
                  DO 220 I = MINI,MAXI
                     IF (IANDJ) MAX = I
                     FACI=FAC*CONI(I)*PNRM(I)*DTWO
                     NN1=NN+1
                     DO 200 J = MINJ,MAX
                        NN = NN+1
                        DIJ(NN)=FACI*CONJ(J)*PNRM(J)
C                    WRITE(6,*) 'WWWDIJ',NN,I,J,II,JJ,DIJ(NN)
  200                CONTINUE
C            CORRECT FOR L-SHELL DOUBLE COUNTING OF THE SP
C            OFF-DIAGONAL TERMS (FOR NON-L SHELLS CSI*CPJ IS ZERO).
C            NN1 POINTS TO THE APPROPRIATE DENSITY ELEMENT
                     IF(MINJ.LE.1.AND.I.GT.1.AND.DOUBLE)
     *                 DIJ(NN1)=DIJ(NN1)*PT5+SPDIJ
  220             CONTINUE
C
C     ----- OVERLAP AND KINETIC ENERGY
C
                  TAA = SQRT(AA1)
                  T1 = -TWO*AJ*AJ*TAA
                  T2 = -PT5*TAA
                  X0 = AX
                  Y0 = AY
                  Z0 = AZ
                  IN = -7
                  DO 320 I = 1,LIT
                     IN = IN+7
                     NI = I
                     DO 300 J = 1,LJT
                        JN = IN+J
                        NJ = J
                        CALL STVINT
                        XIN(JN) = XINT*TAA
                        YIN(JN) = YINT*TAA
                        ZIN(JN) = ZINT*TAA
                        NJ = J+2
                        CALL STVINT
                        XIN(JN+49) = XINT*T1
                        YIN(JN+49) = YINT*T1
                        ZIN(JN+49) = ZINT*T1
                        NJ = J-2
                        IF (NJ .GT. 0) THEN
                           CALL STVINT
                        ELSE
                           XINT = ZERO
                           YINT = ZERO
                           ZINT = ZERO
                        END IF
                        N = (J-1)*(J-2)
                        DUM = N * T2
                        XIN(JN+98) = XINT*DUM
                        YIN(JN+98) = YINT*DUM
                        ZIN(JN+98) = ZINT*DUM
                        IF(LZINT) THEN
                           NJ = J+1
                           CALL STVINT
                           XIN(JN+147) = XINT*TAA
                           YIN(JN+147) = YINT*TAA
C---                       ZIN(JN+147) = ZINT*TAA
                           NJ = J-1
                           IF (NJ .GT. 0) THEN
                              CALL STVINT
                           ELSE
                              XINT = ZERO
                              YINT = ZERO
C---                          ZINT = ZERO
                           END IF
                           XIN(JN+196) = XINT*TAA*NJ
                           YIN(JN+196) = YINT*TAA*NJ
C---                       ZIN(JN+196) = ZINT*TAA*NJ
                        END IF
  300                CONTINUE
  320             CONTINUE
                  DO 340 I = 1,IJ
                     NX = IJX(I)
                     NY = IJY(I)
                     NZ = IJZ(I)
                     YZ = YIN(NY)*ZIN(NZ)
                     DUM = YZ*XIN(NX)
                     DUM1 = (XIN(NX+49)+XIN(NX+98))*YZ+
     *                      (YIN(NY+49)+YIN(NY+98))*XIN(NX)*ZIN(NZ)+
     *                      (ZIN(NZ+49)+ZIN(NZ+98))*XIN(NX)*YIN(NY)
                     SBLK(I) = SBLK(I) + DIJ(I)*DUM
                     TBLK(I) = TBLK(I) + DIJ(I)*(DUM*AJ*FT(I)+DUM1)
                     IF(LZINT) THEN
                        DUM2 = XIN(NX+147)*YIN(NY+196)
     *                       - XIN(NX+196)*YIN(NY+147)
                        ZBLK(I) = ZBLK(I) + DIJ(I)*DUM2*ZIN(NZ)
                     END IF
  340             CONTINUE
C
C     ----- NUCLEAR ATTRACTION
C     EL SIGUIENTE DO SOLO EN CASO NO SCREEN.
C
                  IF (.NOT.SCREEN) THEN
                     DUM = PI212*AA1
                     DO 400 I = 1,IJ
                        DIJ(I) = DIJ(I)*DUM
  400                CONTINUE
                  END IF
C
                  AAX = AA*AX
                  AAY = AA*AY
                  AAZ = AA*AZ
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
                  DO 460 IC = NATST,NATED
                     IF(IC.LE.NAT) THEN
CNEO NEO                --- POSITIVE CHARGE FOR NUCLEAR INTS (NEO) ---
                        ZNUC = ZAN(IC)
                        CX = C(1,IC)
                        CY = C(2,IC)
                        CZ = C(3,IC)
                     ELSE
C
C     COSMO OR CHARMM POINT CHARGES
C
                        IF(ISEPS) THEN
                           ZNUC = -COSZAN(IC-NAT)
                           CX = CORZAN(1,IC-NAT)
                           CY = CORZAN(2,IC-NAT)
                           CZ = CORZAN(3,IC-NAT)
                         ELSE
                           ZNUC = -QCHM(IC-NAT)
                           CX = XCHM(IC-NAT)
                           CY = YCHM(IC-NAT)
                           CZ = ZCHM(IC-NAT)
                        END IF
                     END IF
                     IF(DOESP) THEN
                       CALL FMOATFRG(IC,X(LINDAT),X(LINDATG),X(LIAGLOB),
     *                               X(LIALOC),X(LIABDFG),X(LJABDFG),
     *                               X(LINDBD),X(LFMOZAN),X(LFMOC),
     *                               NATFMO+NBDFG,X(LUNTXYZ),X(LPOPMAT),
     *                               KATOM(II),KATOM(JJ),FRACV,FRACESP,
     *                               KFG,ZNUC,CX,CY,CZ)
C                      WRITE(6,*) 'WWWATOM',IC,FRACV,FRACESP
                       IF(FRACV.EQ.ZERO.AND.FRACESP.EQ.ZERO) GOTO 460
                     ENDIF
C
C         CHECKING IF IT IS AN SCREENED CALCULATION, IF SO CALL SCR1,
C         OTHERWISE, FOLLOW THE GAMESS CODE.
C
                  IF (SCREEN) THEN
                      NN=0
                      DO 425 I=MINI,MAXI
                         L1A=JX(I)
                         M1A=JY(I)
                         N1A=JZ(I)
                         MAX=MAXJ
                         IF(IANDJ) MAX=I
                         DO 415 J=MINJ,MAX
                            NN=NN+1
                            L2B=JX(J)
                            M2B=JY(J)
                            N2B=JZ(J)
                            VAL=FDNAI(VLAMB,AI,AJ,
     *                                L1A,M1A,N1A,L2B,M2B,N2B,
     *                                XI,YI,ZI,XJ,YJ,ZJ,CX,CY,CZ)
                            VBLK(NN) = VBLK(NN) + DIJ(NN)*VAL*ZNUC
C
  415                    CONTINUE
  425                 CONTINUE
                  ELSE
                     XX = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
                     IF (NROOTS.LE.3) CALL RT123
                     IF (NROOTS.EQ.4) CALL ROOT4
                     IF (NROOTS.EQ.5) CALL ROOT5
                     IF (NROOTS.GE.6) CALL ROOT6
                     MM = 0
                     DO 430 K = 1,NROOTS
                        UU = AA*U(K)
                        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 = -7+MM
                        J0 = 2
C
C                       J0 SKIPS THE TRIVIAL COMBINATION NI=NJ=1
C                       THAT IS DONE EXPLICITLY NOW
C
                        XIN(IN+8) = W1
                        YIN(IN+8) = W1
                        ZIN(IN+8) = W1*WW
                        DO 420 I = 1,LIT
                           IN = IN+7
                           NI = I
                           DO 410 J = J0,LJT
                              JN = IN+J
                              NJ = J
                              CALL STVINT
                              XIN(JN) = XINT
                              YIN(JN) = YINT
                              ZIN(JN) = ZINT*WW
  410                      CONTINUE
                           J0 = 1
  420                   CONTINUE
                        MM = MM+49
  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)
                           MM = MM+49
  440                   CONTINUE
                        IF(DOESP) THEN
                           DUMIJ=DUM*DIJ(I)
                           VBLK(I) = VBLK(I) + FRACV*  DUMIJ
                           ZBLK(I) = ZBLK(I) + FRACESP*DUMIJ
                        ELSE
                           VBLK(I) = VBLK(I) + DUM*DIJ(I)
                        ENDIF
  450                CONTINUE
                   END IF
C
  460             CONTINUE
C
C     ----- END OF PRIMITIVE LOOPS -----
C
  500          CONTINUE
  520       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
                  TBLK(NN) = TBLK(NN)*SCAL
                  H(JN) = TBLK(NN) + VBLK(NN)
                  S(JN) = SBLK(NN)
                  T(JN) = TBLK(NN)
                  IF(LZINT) Z(JN) = ZBLK(NN)
                  IF(DOESP) ESP1E(JN)=ZBLK(NN)
  600          CONTINUE
  620       CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
  700    CONTINUE
  720 CONTINUE
C
C     ----- SUM UP PARTIAL CONTRIBUTIONS IF PARALLEL -----
C
      IF (GOPARR) THEN
         CALL DDI_GSUMF(910,H,L2)
         CALL DDI_GSUMF(911,S,L2)
         CALL DDI_GSUMF(912,T,L2)
         IF(LZINT) CALL DDI_GSUMF(913,Z,L2)
C        ONE DDI_GSUMF BELOW CAN BE SPARED
         IF(DOESP) CALL DDI_GSUMF(913,ESP1E,L2)
      END IF
C
C     ----- SAVE H, S, AND T MATRICES ON THE DAF -----
C
      CALL DAWRIT(IDAF,IODA,H,LL2,440,0)
      CALL DAWRIT(IDAF,IODA,S,LL2,441,0)
      CALL DAWRIT(IDAF,IODA,T,LL2,442,0)
C
C     ----- OPTIONAL DEBUG PRINTOUT -----
C
      IF(DBUG) THEN
         WRITE(IW,*) 'OVERLAP MATRIX'
         CALL PRTRIL(S,L1)
         WRITE(IW,*) 'BARE NUCLEUS HAMILTONIAN INTEGRALS (H=T+V)'
         CALL PRTRIL(H,L1)
         WRITE(IW,*) 'KINETIC ENERGY INTEGRALS'
         CALL PRTRIL(T,L1)
         IF(LZINT.OR.DOESP) THEN
            WRITE(IW,*) 'Z-ANGULAR MOMENTUM INTEGRALS'
            CALL PRTRIL(Z,L1)
         END IF
      END IF
      RETURN
C
 9010 FORMAT(1X,'SCREENING VALUE=',F20.5,' IS BEING USED')
      END
C
C*MODULE NEOINT  *DECK ONENI
      SUBROUTINE ONENI
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SOME,DBUG,GOPARR,DSKWRK,MASWRK,LINEAR,LZINT,UNCON,
     *        QUADRE,DOHOP,DOESP1
C
      PARAMETER (MXATM=2000)
      PARAMETER (NDQ=2)
C
      COMMON /FMCOM / X(1)
      COMMON /FMODIM/ MAXBND,MAXKND,MAXCBS,MAXCAO,MAXVEC,MAXL1,
     *                MAXNAT,MAXABD,MAXBAS
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /FMOOPT/ ESPSCA(9),RESPAP(2),RESPPC(2),RESDIM,RESTRI(4),
     *                RCORSD,RESPCT,CONVFG,CNVDMP,COROFF,IXESP,MXITFG,
     *                NGUESS,NBSSE,MODORB,MODPAR,IRSTSTP,IRSTLAY,NPRFMO,
     *                NFMOPAL,MODPRP,MAXL1C,IPIEDA,MODGRD,MODESP,IVMUL
      COMMON /FMOPNT/ LICHFG,LMULFG,LIDMREC,LFRGNAM,LLAYFRG,LINDAT,
     *                LNCBS,LFMOZAN,LFMOC,LFMOMAS,LIZBAS,LIAGLOB,LIBDGH,
     *                LIABDFG,LJABDFG,LNCAO,LIDXCAO,LIAPRJO,LJAPRJO,
     *                LCOREAO,LOCCCOR,LSHIFTB,LIODFMO,LFMODA,LFMODB,
     *                LFMOESPA,LFMOESPB,LLOCFMO,LSCFFRG,LRIJ,LPOPMUL,
     *                LPOPMAT,LIALOC,LINDBD,LIATFRG,LINDFRG,LINDGFRG,
     *                LNATFRG,LNAT0FRG,LIANFRG,LZANFRG,LCFRG,LLIBISH,
     *                LLIBNSH,LLIBNG,LINDATG,LFMOBUF(3),LFMODE,LNUMFRG,
     *                LLOCTAT,LIAOGLOB,LLOADM,LFMOGE,LDGRID,LIODCFMO,
     *                LJOB2GRP,LFMOPG,LEMOCDR,LUNTXYZ,LUNTROT,LSTONEP,
     *                LMAPSU,LFRGMUL
      COMMON /FMORUN/ ESPSCF,E0SCF(2),IDAFMO,ICURLAY,ICURFG,JCURFG,
     *                KCURFG,NCURSH,NGAU,ICURPOP,IFMOSTP,MONCOR,NAT1E,
     *                NEEDR,MODRST,NORBPROJ,NUNESP,ISKIPESP,ICURUNT
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
C
      DATA DEBUG/8HDEBUG   /,DBUGME/8HINT1    /
      DATA ANONE/8HNONE    /
      DATA ANESC/8HNESC    /
C
C     ----- DRIVER FOR ONE ELECTRON INTEGRAL COMPUTATION -----
C     AVOID PRINTING DURING DK NUMERICAL GRADIENT (-2183)
C
      SOME = MASWRK  .AND.  .NOT.(NPRINT.EQ.-5 .OR. NPRINT.EQ.-2183)
      DBUG = MASWRK  .AND.
     *       (NPRINT.EQ.3 .OR. EXETYP.EQ.DEBUG .OR. EXETYP.EQ.DBUGME)
C
      LZINT=LINEAR
      UNCON =RMETHOD.NE.ANONE.AND.MOD(MODQR  ,2).EQ.1
      QUADRE=RMETHOD.NE.ANONE.AND.MOD(MODQR/8,2).EQ.1
      DOHOP=NFG.NE.0.AND.NBDFG.NE.0
C
      IF(SOME) WRITE (IW,9000)
      IF(SOME) CALL TSECND(T0)
C
      L1=NUM
      IF(UNCON) L1=NUMU
      L2 = (NUMNB*NUMNB+NUMNB)/2
C
      CALL VALFM(LOADFM)
      LH   = LOADFM + 1
      LS   = LH     + L2
      LT   = LS     + L2
      LZ   = LT     + L2
      LWRK = LZ     + L2
      IF(.NOT.LZINT) LWRK = LZ
      LAST = LWRK   + L2
      LH0=LAST
      IF(UNCON) LAST=LH0+(NUM*NUM+NUM)/2
      LQQ=LAST
      IF(QUADRE) LAST=LQQ+L2*NDQ*2
      IF(DOHOP) THEN
         LSS=LAST
         LDD=LSS+L1*L1
         LROTLCAO=LDD+(MAXCBS*MAXCBS+MAXCBS)/2
         LQ=LROTLCAO+MAXCBS*MAXCAO
         LSCR=LQ+L1*L1
         LAST=LSCR+L1*8
      ENDIF
      NEED = LAST-LOADFM-1
      CALL GETFM(NEED)
C
C     ----- ORDINARY H, S, AND T INTEGRALS -----
C
      IF(RMETHOD.EQ.ANESC) CALL FLIPBASIS(0)
      IF(UNCON) CALL FLIPBASIS(17)
C
C     QSANDT OVERWRITES THE OVERLAP AND KINETIC ENERGY INTEGRALS
C     (REPLACES THOSE WITH DOUBLE PRECISION BY TRUNCATED QUADRUPLE)
C     PRACTICALLY THERE IS ALMOST NO DIFFERENCE FOR NON-RELATIVISTIC
C     ENERGY BUT LARGE DIFFERENCE FOR RELATIVISTIC (LARGE IN INTEGRALS,
C     THESE ARE MULTIPLIED BY TINY DENSITY SO THE TOTAL ENERGY CHANGES
C     LITTLE, IN THE 12TH DIGIT(?) FOR U).
C
      DOESP1=NFG.NE.0.AND.IFMOSTP.GT.1.AND.IFMOSTP.NE.3.AND.
     *      (NBSSE.NE.3.OR.IFMOSTP.NE.5)
      CALL NUCHST(X(LH),X(LS),X(LT),X(LZ),X(LFMOESPB),
     *            L2,SOME,LZINT,DOESP1,DBUG)
      IF(QUADRE) CALL QSANDT(X(LS),X(LT),X(LQQ),X(LQQ+L2*NDQ),L2,SOME)
      IF(RMETHOD.EQ.ANESC) CALL FLIPBASIS(3)
C
C     ----- DONE WITH NEO 1-PARTICLE INTEGRALS -----
C
      CALL RETFM(NEED)
      IF(MASWRK) WRITE(IW,9090)
      CALL TIMIT(1)
      RETURN
C
 9000 FORMAT(/10X,19(1H*)/10X,'1 NUCLEUS INTEGRALS'/10X,19(1H*))
 9090 FORMAT(1X,'...... END OF ONE-NUCLEUS INTEGRALS ......')
      END
C*MODULE NEOINT  *DECK NJANDK
      SUBROUTINE NJANDK(MIXED,INIT,SCFTYP,DIRNUC,NUCTRN,FA,FB,
     *                  NUCNUC,ELENUC,ELPOVR)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SCHWRZ,PACK2E,PK,PANDK,BLOCK,DIRTRF,FDIFF
      LOGICAL GOPARR,DSKWRK,MASWRK,TDSKWRK,DIRSCF
      LOGICAL INIT,DIRNUC,NUCTRN,MIXED,NUCNUC,ELENUC,ELPOVR
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH)
      PARAMETER (MXAO=8192)
C
      COMMON /FMCOM / XX(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NEOSHL/ NGAUSS,NGAUSE,NGAUSN,NTSHEL,NNUCSH
      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),MIN(MXSH),MAX(MXSH),NSHELL
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PKFIL / PK,PANDK,BLOCK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
C
      DIMENSION FA(*),FB(*)
C     DIMENSION FA(*),FB(*),DA(*),DB(*)
C
      DATA CHECK/8HCHECK   /
C
C     ----- DRIVER FOR CALCULATION OF NUCLEAR INTEGRALS -----
C
C     HONDO INTEGRAL PACKAGE REQUIRES A BUFFER DIMENSIONED FOR
C     THE MAXIMUM ANGULAR MOMENTUM OCCURING IN THE BASIS SET.
C     A PURE SP BASIS USING THE POPLE INTEGRALS DOESN'T NEED THIS.
C
      NSHLSV = NSHELL
      NSHELL = NTSHEL
      CALL BASCHK(LMAX)
      NSHELL = NSHLSV
                    NANGM =  4
      IF(LMAX.EQ.2) NANGM =  6
      IF(LMAX.EQ.3) NANGM = 10
      IF(LMAX.EQ.4) NANGM = 15
                 MAXG = NANGM**4
      IF(PK)     MAXG = MAXG*3
C
C         NEXT IS TO MAKE PLACEHOLDER BETA STORAGE APPEAR TO BE USED
C
      FB(1)=0.0D+00
C
      IF(MIXED) THEN
         NSH2 = NSHELL*NNUCSH
      ELSE
         NSH2 = (NNUCSH*NNUCSH+NNUCSH)/2
      END IF
C
      CALL VALFM(LOADFM)
C
      IF(DIRNUC) THEN
         LXINTS = LOADFM + 1
         LGHOND = LXINTS + NSH2
         LDDIJ  = LGHOND + MAXG
         LDSH   = LDDIJ  + 49*MXG2
         LAST   = LDSH   + NSH2
         LBUFP  = LAST
         LBUFK  = LAST
         LIX    = LAST
      ELSE
         LBUFP  = LOADFM + 1
         LBUFK  = LBUFP  + NINTMX
         LIX    = LBUFK  + NINTMX
         LXINTS = LIX    + NINTMX
         LGHOND = LXINTS + NSH2
         LDDIJ  = LGHOND + MAXG
         LDSH   = LDDIJ  + 49*MXG2
         LAST   = LDSH
      END IF
C
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      TDSKWRK = DSKWRK
      DSKWRK  = .TRUE.
C
C     ---- INITIALIZATION FOR ELECTRON-NUCLEUS OR NUCLEUS-NUCLEUS
C          INTEGRALS ----
C
      IF(INIT) THEN
         CALL NDEBUT(DIRSCF,XX(LBUFP),XX(LBUFK),XX(LIX),
     *               NINTMX,NEED,DIRTRF,FDIFF,MIXED)
         IF(EXETYP.EQ.CHECK) GOTO 200
         IF(DIRNUC.OR.NUCTRN) GOTO 200
      END IF
C
C     ---- CALCULATE INTEGRALS REQUIRED FOR THE SCHWARTZ INEQUALITY ----
C     ******** PRESENTLY SCHWARTZ INEQUALITY NOT USED FOR NEO INTEGRALS
C
C     ---- NOW CALCULATE THE INTEGRALS - FOR DIRECT OPTION ----
C          CONTRIBUTIONS TO THE FOCK MATRICES ARE ADDED ON-THE-FLY
C
      SCHWRZ = .FALSE.
      CALL TIMIT(1)
      N1=1
      N2=1
      DUM=0.0D+00
      IF(ELPOVR) THEN
         IDFTOVRLP=4
         CALL NEODFTINT(IDFTOVRLP)
      ELSE
         CALL TWOENI(SCFTYP,DIRNUC,NUCTRN,INTTYP,SCHWRZ,NINT,NSCHWZ,
     *           N1,N2,XX(LBUFP),XX(LBUFK),XX(LIX),NINTMX,XX(LXINTS),
     *           NSH2,XX(LGHOND),MAXG,XX(LDDIJ),IA,DUM,FA,
     *           DUM,NTSHEL,MIXED,NUCNUC,ELENUC)
      END IF
C
  200 CONTINUE
C
      DSKWRK  = TDSKWRK
C
      CALL RETFM(NEED)
C
      IF (MASWRK .AND. INIT) THEN
         IF(MIXED.AND.ELENUC.AND..NOT.ELPOVR) THEN
           WRITE(IW,*)' .. END OF ELECTRON-NUCLEAR COULOMB INTEGRALS ..'
         ELSE IF(MIXED.AND..NOT.ELENUC.AND.ELPOVR) THEN
         WRITE(IW,*)' '
         ELSE
           WRITE(IW,*)' ....... END OF NUCLEAR-NUCLEAR INTEGRALS ......'
         END IF
      END IF
      CALL TEXIT(1,1)
C
      RETURN
      END
C
C*MODULE NEOINT  *DECK TWOENI
      SUBROUTINE TWOENI(SCFTYP,DIRSCF,DIRTRF,
     *                 INTTYP,SCHWRZ,NINT,NSCHWZ,L1,L2,
     *                 BUFP,BUFK,IX,NINTMX,
     *                 XINTS,NSH2,GHONDO,MAXG,DDIJ,
     *                 IA,DA,FA,DSH,
     *                 NTSHEL,MIXED,NUCNUC,ELENUC)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DIRSCF,DIRTRF,CMBDIR
      LOGICAL OUT,SCHWRZ,SCHSKP,GOPARR,DSKWRK,MASWRK,NXT
      LOGICAL SKIPA,SKIPB,SKIPC,NPSYM
      LOGICAL PK,PANDK,NOTPK,BLOCK,POPLE,GPSAVE
      LOGICAL NUCNUC,ELENUC,MIXED,FIRST
C
      LOGICAL DIRNUC,NUFOCK
C
      DIMENSION BUFP(NINTMX),BUFK(NINTMX),IX(*),XINTS(NSH2),
     *          GHONDO(MAXG),IA(L1),DA(L2),FA(L2),
     *          DSH(NSH2),DDIJ(*)
C
      DIMENSION MI(48),MJ(48),MK(48),M0(48),NORGP(3),NORGH(3)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
      COMMON /GOUT  / GPOPLE(768),NORG
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
C
      COMMON /NUCMOP/ DIRNUC,NUFOCK
C
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PKFIL / PK,PANDK,BLOCK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /SHLG70/ ISH,JSH,KSH,LSH,IJKLXX(4)
      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
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
      PARAMETER (ZERO=0.0D+00)
C
      DATA NORGP/0,  256,   512/
C
C     MANY ARGUMENTS ARE OPTIONAL, YOU MUST ALLOCATE STORAGE FOR
C     ALL CALLS            : GHONDO, DDIJ, XINTS
C     CONVENTIONAL         : BUFP, IX, AND POSSIBLY BUFK
C     DIRECT SCF           : IA, DSH, DA, FA
C     DIRECT NLO           : DNLO, FNLO
C     DIRECT TRANSFORMATION: BUFP, IX
C
C     ----- NUCLEAR-NUCLEAR INTEGRALS OR       -----
C     ----- ELECTRON-NUCLEAR COULOMB INTEGRALS -----
C     ----- THIS VERSION CAN HANDLE G SHELLS   -----
C
      TIM = ZERO
      CALL TSECND(TIM)
C
      FIRST=.TRUE.
C
C     ----- INITIALIZATION FOR PARALLEL WORK -----
C
      NXT = IBTYP.EQ.1
      IPCOUNT = ME - 1
      NEXT = -1
      MINE = -1
      ICONT=  0
      CMBDIR= DIRSCF .OR. DIRTRF
      GPSAVE = GOPARR
      IF(DIRTRF) GOPARR=.FALSE.
C
      NSHLSV = NSHELL
      NSHELL = NTSHEL
      CALL BASCHK(LMAX)
      NSHELL = NSHLSV
C
                    NANGM =  4
      IF(LMAX.EQ.2) NANGM =  6
      IF(LMAX.EQ.3) NANGM = 10
      IF(LMAX.EQ.4) NANGM = 15
      NORGH(1) = 0
      NORGH(2) = NORGH(1) + NANGM**4
      NORGH(3) = NORGH(2) + NANGM**4
C
      NOTPK = .NOT.PK
      NINT  = 0
      NSCHWZ= 0
      SCHSKP=.FALSE.
      DENMAX = ZERO
C
C     ---- SET LOOP LIMITS FOR TYPE OF INTEGRAL REQUIRED ----
C
      I0 = NSHELL + 1
      J0 = NSHELL + 1
      IF(MIXED) THEN
         K0 = 1
         L0 = 1
      ELSE
         K0 = NSHELL + 1
         L0 = NSHELL + 1
      END IF
      II = I0
      JJ = J0
      KK = K0
      LL = L0
C
C     ----- I SHELL -----
C
      DO 920 II = I0,NTSHEL
C
C     ----- CHECK CPU TIME -----
C
      CALL TSECND(TIM)
      IF(TIM.GE.TIMLIM) THEN
         IF(MASWRK) THEN
            IF(MIXED) THEN
               WRITE(IW,9030)
             ELSE
               WRITE(IW,9040)
             END IF
         END IF
         WRITE(IW,9050) SCFTYP
         CALL ABRT
         STOP
      END IF
C
C     ----- PRINT INTERMEDIATE RESTART DATA -----
C
      IF(NPRINT.NE.-5  .AND.  .NOT.CMBDIR .AND. MASWRK)
     *   WRITE(IW,9010) II,JST,KST,LST,NREC,ICOUNT
C
C     ----- SKIP I SHELL IF NOT SYMMETRY UNIQUE -----
C
      DO 120 IT = 1,NT
         ID = MAPSHL(II,IT)
         IF (ID .GT. II) GO TO 920
         MI(IT) = ID
  120 CONTINUE
C
C     ---- OBTAIN SMALLEST EXPONENT IN THIS CONTRACTION -----
C
      ZMINI = 1.0D+30
      IG1 = KSTART(II)
      IG2 = IG1 + KNG(II) - 1
      DO IG = IG1,IG2
         ZMINI =MIN(ZMINI,EX(IG))
      END DO
C
C     ----- J SHELL -----
C
      DO 900 JJ = J0,II
         DO 200 IT = 1,NT
            ID = MI(IT)
            JD = MAPSHL(JJ,IT)
            MJ(IT) = JD
            IF (ID .GE. JD) GO TO 160
            ND = ID
            ID = JD
            JD = ND
  160       IF (ID-II) 200,180,900
  180       IF (JD-JJ) 200,200,900
  200    CONTINUE
C
C        ----- GO PARALLEL! -----
C
         IF (NXT .AND. GOPARR) THEN
            MINE = MINE + 1
            IF (MINE.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
            IF (NEXT.NE.MINE) GO TO 900
         END IF
C
         ZMINJ = 1.0D+30
         JG1 = KSTART(JJ)
         JG2 = JG1 + KNG(JJ) - 1
         DO JG = JG1,JG2
            ZMINJ = MIN(ZMINJ,EX(JG))
         END DO
C
C        ----- K SHELL -----
C
         IF(MIXED) THEN
            LIMK = NSHELL
         ELSE
            LIMK = JJ
         END IF
         DO 880 KK = K0,LIMK
C           IATOM = KATOM(KK)
C           IF(IATOM.EQ.2) GOTO 880
            DO 340 IT = 1,NT
               ID = MI(IT)
               JD = MJ(IT)
               KD = MAPSHL(KK,IT)
               MK(IT) = KD
  240          IF (ID .GE. JD) GO TO 260
               ND = ID
               ID = JD
               JD = ND
  260          CONTINUE
C              IF(.NOT.MIXED) THEN
                  IF (JD .GE. KD) GO TO 280
                  ND = JD
                  JD = KD
                  KD = ND
                  GO TO 240
C              END IF
  280          IF (ID-II) 340,300,880
  300          IF (JD-JJ) 340,320,880
  320          IF (KD-KK) 340,340,880
  340       CONTINUE
C
            ZMINK = 1.0D+30
            KG1 = KSTART(KK)
            KG2 = KG1 + KNG(KK) - 1
            DO KG = KG1,KG2
               ZMINK = MIN(ZMINK,EX(KG))
            END DO
C
C           ----- L SHELL ----
C
            DO 860 LL = L0,KK
               N4 = 0
               DO 540 IT = 1,NT
               ID = MI(IT)
               JD = MJ(IT)
               KD = MK(IT)
               LD = MAPSHL(LL,IT)
  380          IF (ID .GE. JD) GO TO 400
               ND = ID
               ID = JD
               JD = ND
  400          CONTINUE
C              IF(.NOT.MIXED) THEN
                  IF (JD .GE. KD) GO TO 420
                  ND = JD
                  JD = KD
                  KD = ND
                  GO TO 380
C              END IF
  420          IF (KD .GE. LD) GO TO 440
               ND = KD
               KD = LD
               LD = ND
               GO TO 400
  440          IF (ID-II) 540,460,860
  460          IF (JD-JJ) 540,480,860
  480          IF (KD-KK) 540,500,860
  500          IF (LD-LL) 540,520,860
  520          N4 = N4+1
               M0(N4) = IT
  540          CONTINUE
C
C              ----- GO PARALLEL! -----
C
               IF ((.NOT.NXT) .AND. GOPARR) THEN
                  IPCOUNT = IPCOUNT + 1
                  IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 860
               END IF
C
               ZMINL=1.0D+30
               LG1=KSTART(LL)
               LG2=LG1+KNG(LL)-1
               DO LG=LG1,LG2
                  ZMINL=MIN(ZMINL,EX(LG))
               END DO
               ZMIN = MIN(ZMINI,ZMINJ,ZMINK,ZMINL)
C
C              ----- CHECK FOR REDUNDANIES BETWEEN THE 3 COMBINATIONS
C                    (II,JJ//KK,LL),(II,KK//JJ,LL),(II,LL//JJ,KK)
C
               IF(MIXED) THEN
                  SKIPA = .FALSE.
                  SKIPB = .FALSE.
               ELSE
                  SKIPA =  JJ.EQ.KK
                  SKIPB = (II.EQ.KK) .OR. (JJ.EQ.LL)
               END IF
               SKIPC = (II.EQ.JJ) .OR. (KK.EQ.LL)
               NPSYM = .FALSE.
               IF (SKIPA .OR. SKIPB .OR. SKIPC) GO TO 720
               NPSYM = .TRUE.
               DO 640 M = 1,N4
                  IT = M0(M)
                  IH = MI(IT)
                  JH = MJ(IT)
                  IF(JH.LE.IH) THEN
                     ID = IH
                     JD = JH
                  ELSE
                     ID = JH
                     JD = IH
                  END IF
                  IF(.NOT.MIXED) THEN
                     IF(.NOT.SKIPA) SKIPA =
     *                  (ID.EQ.II .AND. JD.EQ.KK) .OR.
     *                  (ID.EQ.JJ .AND. JD.EQ.LL)
                     IF(.NOT.SKIPB) SKIPB =
     *                  (ID.EQ.II .AND. JD.EQ.LL) .OR.
     *                  (ID.EQ.JJ .AND. JD.EQ.KK)
                  END IF
                  IF (SKIPA .AND. SKIPB) GO TO 660
                  KH = MK(IT)
                  IF(.NOT.MIXED) THEN
                  IF(KH.LE.IH) THEN
                     ID = IH
                     KD = KH
                  ELSE
                     ID = KH
                     KD = IH
                  END IF
                  IF(.NOT.SKIPC) SKIPC = (ID.EQ.II .AND. KD.EQ.LL) .OR.
     *                                   (ID.EQ.JJ .AND. KD.EQ.KK)
                  END IF
                  IF(SKIPA .AND. SKIPC) GO TO 680
                  IF(SKIPB .AND. SKIPC) GO TO 700
  640             CONTINUE
                  GO TO 720
C
  660             SKIPC = .TRUE.
                  GO TO 720
  680             SKIPB = .TRUE.
                  GO TO 720
  700             SKIPA = .TRUE.
C
  720             CONTINUE
                  Q4 = NT
                  Q4 = Q4 / N4
C
C                 ----- (II,JJ//KK,LL) -----
C
                  IEXCH = 1
                  ISH = II
                  JSH = JJ
                  KSH = KK
                  LSH = LL
                  QQ4 = Q4
                  IF(SKIPA .AND. NPSYM) QQ4 = QQ4+Q4
                  IF(SKIPB .AND. NPSYM) QQ4 = QQ4+Q4
                  GO TO 780
C
C                 ----- (II,KK//JJ,LL) -----
C
  740             IF (SKIPA) GO TO 760
                  IF (MIXED) GO TO 840
                  IEXCH = 2
                  ISH = II
                  JSH = KK
                  KSH = JJ
                  LSH = LL
                  QQ4 = Q4
                  IF (SKIPC .AND. NPSYM) QQ4 = QQ4+Q4
                  GO TO 780
C
C                 ----- (II,LL//JJ,KK) -----
C
  760             IF (SKIPB .OR. SKIPC) GO TO 840
                  IF (MIXED) GO TO 840
                  IEXCH = 3
                  ISH = II
                  JSH = LL
                  KSH = JJ
                  LSH = KK
                  QQ4 = Q4
C
C                 ----- COMPUTE TWO-ELECTRON INTEGRALS ----
C
  780             CONTINUE
C
C                 APPLY THE SCHWARZ INEQUALITY, WHICH IS
C                 (II,JJ//KK,LL).LE. SQRT((II,JJ//II,JJ)*(KK,LL//KK,LL))
C                 SEE J.L.WHITTEN, J.CHEM.PHYS. 58,4496-4501(1973)
C
                  IF(SCHWRZ) THEN
                     IJIJ = (ISH*ISH-ISH)/2 + JSH
                     KLKL = (KSH*KSH-KSH)/2 + LSH
                     TEST = QQ4*XINTS(IJIJ)*XINTS(KLKL)
                     IF(DIRSCF) THEN
                        DENMAX = SCHWDN(DSH,ISH,JSH,KSH,LSH,IA)
                        TEST = TEST*DENMAX
                     END IF
                     SCHSKP = TEST.LT.CUTOFF
                     IF(SCHSKP) NSCHWZ = NSCHWZ + 1
                  END IF
C
C                 USE POPLE CODE FOR ANY PURE SP INTEGRAL BLOCKS,
C                 USE HONDO RYS POLYNOMIAL CODE FOR OTHER BLOCKS
C
C      MWS - NOTE THAT INTTYP=0 CORRESPONDS TO BEST TIMING,
C            1,2,3=VARIOUS SUBSETS (SINCE 11/2004).
C      SO, THE 2ND LINE HERE FORCES RYS QUADRATURE BE USED, ALWAYS.
C
C      THE GENR70 CALL PROBABLY NEEDS CONSTANT ARGUMENTS TO WORK,
C      THE ZMIN STUFF CAN BE DELETED (ROTATED AXIS NOW ACCURATE),
C      AND SO ON (COMPARE TO TWOEI AND SHELLQUART).
C
                  POPLE=.TRUE.
                  IF(INTTYP.GE.0) POPLE=.FALSE.
                  IF(KTYPE(ISH).GT.2) POPLE=.FALSE.
                  IF(KTYPE(JSH).GT.2) POPLE=.FALSE.
                  IF(KTYPE(KSH).GT.2) POPLE=.FALSE.
                  IF(KTYPE(LSH).GT.2) POPLE=.FALSE.
                  IF(ZMIN.LT.0.10D+00) POPLE=.FALSE.
C
                  IF(POPLE) THEN
                     NORG = NORGP(IEXCH)
                     CALL GENR70(IEXCH,SCHSKP)
                  ELSE
                     IF(PK  .AND.  IEXCH.EQ.1)
     *               CALL ZPKOUT(ISH,JSH,KSH,LSH,GHONDO,
     *                           SKIPA,SKIPB,SKIPC,NPSYM)
                     IF(SCHSKP) GO TO 820
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
                     NSHLSV = NSHELL
                     NSHELL = NTSHEL
                     NORG=0
                     IF(PK) NORG=NORGH(IEXCH)
                     CALL NEO_SHELLS(1,ISH,JSH,KSH,LSH,FIRST)
                     CALL IJPRIM(DDIJ)
                     IF(NIJ.EQ.0) THEN
                        NSHELL = NSHLSV
                        GO TO 820
                     END IF
                     CALL NEO_SHELLS(2,ISH,JSH,KSH,LSH,FIRST)
                     IF(NOTPK) CALL NEO_ZQOUT(GHONDO,MIXED)
                     NSHELL = NSHLSV
C
C               ----- DO INTEGRAL BATCH, SSSS IS A SPECIAL CASE -----
C
                     IF(IJKL.EQ.1) THEN
                        CALL S0000(GHONDO,DDIJ)
                     ELSE
                        CALL GENRAL(GHONDO,DDIJ)
                     END IF
                  END IF
C
                  IF(DIRNUC) THEN
                     IF(NUCNUC) THEN
                        CALL NNFOCK(IEXCH,POPLE,IA,DA,FA,GHONDO,NINT)
                     ELSE
                        CALL DDFOCK(IEXCH,POPLE,IA,DA,FA,GHONDO,NINT)
                     ENDIF
C
                  ELSE IF(DIRTRF) THEN
C                        NEXT WON'T WORK IF POPLE INTEGRALS GENERATED
                     CALL DIRTRN(BUFP,IX,NINTMX,GHONDO,ICONT,NINT)
                  ELSE
                     IF(NOTPK) CALL NEO_QOUT(IEXCH,POPLE,BUFP,IX,NINTMX,
     *                                   GHONDO,MIXED)
                  END IF
C
  820             CONTINUE
                  GO TO (740,760,840),IEXCH
C
C                 ----- WRITE THE PK INTEGRALS TO FILE (IS) -----
C
  840             CONTINUE
C                        NEXT WON'T WORK IF POPLE INTEGRALS GENERATED
                  IF(PK) CALL PKFILE(II,JJ,KK,LL,SKIPA,SKIPB,SKIPC,
     *                               NPSYM,BUFP,BUFK,IX,NINTMX,GHONDO)
  860          CONTINUE
  880       CONTINUE
  900    CONTINUE
  920 CONTINUE
C
C     ----- END OF SHELL LOOPS -----
C
      IF(NXT  .AND.  GOPARR) CALL DDI_DLBRESET
C
C     ----- OUTPUT THE LAST BITS OF INTEGRALS -----
C
      IF(DIRTRF) CALL ONEIDX(BUFP,IX,ICONT)
      IF(.NOT.CMBDIR) THEN
         IF(SCHWRZ) THEN
            IF(GOPARR) CALL DDI_GSUMI(1055,NSCHWZ,1)
            IF(NPRINT.NE.-5 .AND. MASWRK) WRITE(IW,9020) NSCHWZ
         END IF
         CALL NEO_FINAL(1,II,II,II,II,PANDK,BUFP,BUFK,IX,NINTMX,
     *              NUCNUC,ELENUC,MIXED)
      END IF
      GOPARR = GPSAVE
      RETURN
C
 9010 FORMAT(1X,'II,JST,KST,LST =',4I3,' NREC =',I10,' INTLOC =',I5)
 9020 FORMAT(1X,'SCHWARZ INEQUALITY TEST SKIPPED',I10,
     *        ' INTEGRAL BLOCKS.')
 9030 FORMAT(//1X,'*** THIS JOB HAS EXHAUSTED ITS CPU TIME ***'/
     *         1X,'     (WHILE COMPUTING E-P+ INTEGRALS)'///)
 9040 FORMAT(//1X,'*** THIS JOB HAS EXHAUSTED ITS CPU TIME ***'/
     *         1X,'     (WHILE COMPUTING 2P+ INTEGRALS)'///)
 9050 FORMAT(1X,'PLEASE RESTART THIS NEO-',A8,' CALCULATION')
      END
C
C*MODULE NEOINT  *DECK NDEBUT
      SUBROUTINE NDEBUT(DIRSCF,BUFP,BUFK,IX,NINTMX,NEED,
     *                 DIRTRF,FDIFF,ELNU)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL PK,PANDK,BLOCK,OUT,DIRSCF,DIRTRF,FDIFF,DIR
      LOGICAL GOPARR,DSKWRK,MASWRK,ELNU
C
      DIMENSION BUFP(NINTMX),BUFK(NINTMX),IX(*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000,  MXATM=2000)
      PARAMETER (MXAO=8192)
C
      COMMON /B     / X(MXSH),Y(MXSH),Z(MXSH)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTPR / QINT(2),VALINT(2),JCINT(11)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /NEOSHL/ NGAUSS,NGAUSE,NGAUSN,NTSHEL,NNUCSH
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
      COMMON /PICON / PITO52,PIDIV4
      COMMON /PKFIL / PK,PANDK,BLOCK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
C
      PARAMETER (TEN=1.0D+01, ONE=1.0D+00, TWO=2.0D+00,
     *           TWOPT5=2.5D+00, FOUR=4.0D+00, E=2.30258D+00)
C
      DATA CHECK/8HCHECK   /
C
C     ----- INITIALIZE TWO ELECTRON INTEGRAL CALCULATION -----
C
      DIR = DIRSCF .OR. DIRTRF
C
      PIDIV4 = ATAN(ONE)
      PI = FOUR*PIDIV4
      PITO52 = TWO * PI**TWOPT5
C
      TOL = E*ITOL
      CUTOFF = ONE/(TEN**ICUT)
      IF(FDIFF) CUTOFF = CUTOFF/TWO
C
      DO 100 I=1,MXAO
         IA(I) = (I*I-I)/2
  100 CONTINUE
C
C     ----- ERIC (AND OTHER INTEGRAL PACKAGE) INITIALIZATIONS -----
C
      CALL ERIPRE
C
      IF((NPRINT.NE.-5) .AND. MASWRK) THEN
         IF(DIR) THEN
C
            IF(ELNU) THEN
               WRITE(IW,9000)
            ELSE
               WRITE(IW,9001)
            ENDIF
C
            IF(DIRSCF) WRITE(IW,9040)
            IF(DIRTRF) WRITE(IW,9045)
         ELSE
            IF(ELNU) THEN
               WRITE(IW,9000)
            ELSE
               WRITE(IW,9001)
            END IF
            IF(PK) THEN
               IF(     PANDK) WRITE(IW,9010)
               IF(.NOT.PANDK) WRITE(IW,9020)
            ELSE
               WRITE(IW,9030)
            END IF
            NBYTES = 8 + (LABSIZ*8)/NWDVAR
            IF(PANDK) NBYTES = NBYTES+8
            WRITE(IW,9050) NINTMX,NBYTES
            IF (ELNU) THEN
               WRITE(IW,9061) NEED
            ELSE
               WRITE(IW,9060) NEED
            END IF
         END IF
      END IF
C
      OUT = NPRINT.EQ.4 .AND. MASWRK
      JCINT(1) = 0
C
      DO 200 I = 1,NTSHEL
         ICC = KATOM(I)
         X(I) = C(1,ICC)
         Y(I) = C(2,ICC)
         Z(I) = C(3,ICC)
  200 CONTINUE
C
      IF(EXETYP.EQ.CHECK)  GO TO 400
      IF(DIRSCF.OR.DIRTRF) GO TO 400
C
      CALL SEQREW(IS)
C
      IF(IREST.LT.1) GO TO 400
      IF(NREC.LE.1  .OR.  INTLOC.LE.1) GO TO 400
C
C     ----- POSITION THE INTEGRAL FILE FOR A RESTART JOB -----
C
      ICOUNT = INTLOC
      N = NREC-1
      IF (MASWRK) THEN
         DO 310 I = 1,N
            READ(IS)
  310    CONTINUE
      END IF
      IF(.NOT.PANDK) CALL PREAD (IS,BUFP,IX,NXX,NINTMX)
      IF(     PANDK) CALL PKREAD(IS,BUFP,BUFK,IX,NXX,NINTMX)
      CALL SEQREW(IS)
      IF (MASWRK) THEN
         DO 320 I = 1,N
            READ(IS)
  320    CONTINUE
      END IF
      RETURN
C
C     ----- NORMAL START -----
C
  400 CONTINUE
      IF(IST.LE.1) IST = 1
      IF(JST.LE.1) JST = 1
      IF(KST.LE.1) KST = 1
      IF(LST.LE.1) LST = 1
      NREC   = 1
      INTLOC = 1
      ICOUNT = 1
      RETURN
C
 9000 FORMAT(/10X,26(1H-)/10X,'ELECTRON-NUCLEUS INTEGRALS'/10X,26(1H-)/)
 9001 FORMAT(/10X,19(1H-)/10X,'2 NUCLEUS INTEGRALS'/10X,19(1H-)/)
 9010 FORMAT(1X,'THE -PK- OPTION IS ON, CREATING -P- AND -K-',
     *          ' SUPERMATRICES.')
 9020 FORMAT(1X,'THE -PK- OPTION IS ON, CREATING A -P- SUPERMATRIX.')
 9030 FORMAT(1X,'THE -PK- OPTION IS OFF, THE INTEGRALS ARE NOT IN',
     *          ' SUPERMATRIX FORM.')
 9040 FORMAT(1X,'DIRECT SCF METHOD SKIPS INTEGRAL STORAGE ON DISK.')
 9045 FORMAT(1X,'DIRECT TRANSFORMATION SKIPS AO INTEGRAL STORAGE',
     *          ' ON DISK.')
 9050 FORMAT(1X,'STORING',I8,' INTEGRALS/RECORD ON DISK, USING',I3,
     *          ' BYTES/INTEGRAL.')
 9060 FORMAT(1X,'NUCLEAR-NUCLEAR INTEGRAL EVALUATION REQUIRES',I8,
     *          ' WORDS OF MEMORY.')
 9061 FORMAT(1X,'ELECTRON-NUCLEAR INTEGRAL EVALUATION REQUIRES',I8,
     *          ' WORDS OF MEMORY.')
      END
C
C*MODULE NEOINT  *DECK NEO_FINAL
      SUBROUTINE NEO_FINAL(INDEX,II,JJ,KK,LL,PANDK,BUFP,BUFK,
     *                 IX,NINTMX,NUCNUC,ELENUC,MIXED)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL OUT,PANDK,GOPARR,DSKWRK,MASWRK,NUCNUC,ELENUC,MIXED
C
      DIMENSION BUFP(NINTMX),BUFK(NINTMX),IX(*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      COMMON /INTPR / Q(2),V(2),JC,N1(2),J1(2),J2(2),J3(2),J4(2)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /SHLT  / TOL,CUTOFF,ICOUNT,OUT
C
      IF(OUT .AND. JC.GT.0)
     * WRITE(IW,9000) (J1(M),J2(M),J3(M),J4(M),Q(M),N1(M),V(M),M=1,JC)
C
      IF (INDEX.EQ.1) GO TO 140
C
C         RAN OUT OF TIME, FLUSH PARTIAL BUFFER, PRINT RESTART DATA
C
      IREST = 1
      IST = II
      JST = JJ
      KST = KK
      LST = LL+1
      IF(LST.LE.KK) GO TO 120
      LST = 1
      KST = KK+1
      IF(KST.LE.JJ) GO TO 120
      KST = 1
      JST = JJ+1
      IF(JST.LE.II) GO TO 120
      JST = 1
      IST = II+1
      IF(IST.GT.NSHELL) GO TO 140
C
  120 CONTINUE
      NXX = ICOUNT-1
      IF(.NOT.PANDK) CALL PWRIT (IS,BUFP,IX,NXX,NINTMX)
      IF(     PANDK) CALL PKWRIT(IS,BUFP,BUFK,IX,NXX,NINTMX)
      NINT = NINTMX*(NREC-1)+ICOUNT-1
      IF (MASWRK) THEN
         WRITE(IW,9010) NINT,NREC,IS
         WRITE(IW,9020) NREC,ICOUNT,IST,JST,KST,LST
      END IF
      RETURN
C
C        DONE WITH INTEGRALS, WRITE LAST BUFFER, PRINT STATISTICS
C
  140 CONTINUE
      IREST = 0
      IST = 1
      JST = 1
      KST = 1
      LST = 1
      NXX = ICOUNT-1
      NXX = -NXX
      IF(.NOT.PANDK) CALL PWRIT (IS,BUFP,IX,NXX,NINTMX)
      IF(     PANDK) CALL PKWRIT(IS,BUFP,BUFK,IX,NXX,NINTMX)
      NINT = NINTMX*(NREC-1)+ICOUNT-1
      IF (GOPARR) THEN
         CALL DDI_GSUMI(1056,NINT,1)
         CALL DDI_GSUMI(1057,NREC,1)
      END IF
      IF (MASWRK) THEN
         IF(ELENUC .AND. MIXED) THEN
            WRITE(IW,9010) NINT,NREC,IS
         ELSE IF (NUCNUC) THEN
            WRITE(IW,9011) NINT,NREC,IS
         ELSE
            WRITE(IW,9012) NINT,NREC,IS
         END IF
      END IF
      RETURN
C
 9000 FORMAT(3(4I3,F6.3,I5,E20.12))
 9010 FORMAT(1X,'TOTAL NUMBER OF NONZERO ELECTRON-NUCLEAR INTEGRALS =',
     *       I20/
     *       1X,I10,' INTEGRAL RECORDS WERE STORED ON DISK FILE',I3,'.')
 9011 FORMAT(1X,'TOTAL NUMBER OF NONZERO NUCLEAR-NUCLEAR INTEGRALS =',
     *       I20/
     *       1X,I10,' INTEGRAL RECORDS WERE STORED ON DISK FILE',I3,'.')
 9012 FORMAT(1X,'TOTAL NUMBER OF NONZERO TWO-ELECTRON INTEGRALS =',I20/
     *       1X,I10,' INTEGRAL RECORDS WERE STORED ON DISK FILE',I3,'.')
 9020 FORMAT(/1X,'...... WARNING   .......   WARNING   .......'/
     *        1X,'TIME LIMIT HAS EXPIRED.  THIS JOB MUST BE RESTARTED.'/
     *        1X,'IF YOU SAVED THE INTEGRALS, RESTART WITH IREST=1,'/
     *        1X,'NREC=',I8,' INTLOC=',I6,' IST,JST,KST,LST=',4I6)
      END
C
C*MODULE NEOINT  *DECK NEO_QOUT
      SUBROUTINE NEO_QOUT(IEXCH,POPLE,BUFP,IX,NINTMX,GHONDO,MIXED)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL IANDJ,KANDL,SAME,OUT,POPLE,MIXED
C
      DIMENSION BUFP(NINTMX),IX(*),GHONDO(*)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      COMMON /FLIPS / IB(4,3)
      COMMON /GOUT  / GPOPLE(768),NORG
      COMMON /INTDEX/ IJGT(784),IJX(784),IJY(784),IJZ(784),IK(784),
     *                KLGT(784),KLX(784),KLY(784),KLZ(784)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MISC  / IANDJ,KANDL,SAME
      COMMON /NSHEL / EXX(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 /PCKLAB/ LABSIZ
      COMMON /SHLG70/ ISH,JSH,KSH,LSH,IJKLXX(4)
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      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 IBPOP(4,4)
C
C-NEXT STATEMENT IS FOR VARIOUS IBM XLF 3.X AND 5.X COMPILERS-
C
      SAVE IJN,KLN,IBB,JBB,KBB,LBB
C
      DATA IBPOP/0,0,0,0,64,16,4,1,128,32,8,2,192,48,12,3/
      DATA HALF /0.5D+00/
C
C     ----- PACK THE 4 INDICES OF INTEGRAL INTO ONE WORD
C     ----- WRITE LABEL + INTEGRAL ON TAPE (IS)
C
      IF(POPLE) THEN
         IANDJ = ISH.EQ.JSH
         KANDL = KSH.EQ.LSH
         SAME  = ISH.EQ.KSH  .AND.  JSH.EQ.LSH
         IBB  = IB(1,IEXCH)
         JBB  = IB(2,IEXCH)
         KBB  = IB(3,IEXCH)
         LBB  = IB(4,IEXCH)
         MINI = KMIN(ISH)
         MINJ = KMIN(JSH)
         MINK = KMIN(KSH)
         MINL = KMIN(LSH)
         MAXI = KMAX(ISH)
         MAXJ = KMAX(JSH)
         MAXK = KMAX(KSH)
         MAXL = KMAX(LSH)
         LOCI = KLOC(ISH)-MINI
         LOCJ = KLOC(JSH)-MINJ
         LOCK = KLOC(KSH)-MINK
         LOCL = KLOC(LSH)-MINL
      END IF
C
      IJN = 0
      JMAX = MAXJ
      DO 260 I = MINI,MAXI
         IF (IANDJ) JMAX = I
         DO 240 J = MINJ,JMAX
            IJN = IJN+1
            IF(POPLE) THEN
               N1 = IBPOP(IBB,I)+IBPOP(JBB,J)+1+NORG
            ELSE
               N1 = IJGT(IJN)
            END IF
            LMAX = MAXL
            KLN = 0
            DO 220 K =  MINK,MAXK
               IF (KANDL) LMAX = K
               DO 200 L = MINL,LMAX
                  KLN = KLN+1
                  IF(.NOT.MIXED) THEN
                     IF(SAME  .AND.  KLN.GT.IJN) GO TO 240
                  END IF
                  IF(POPLE) THEN
                     NN = N1+IBPOP(KBB,K)+IBPOP(LBB,L)
                     VAL = GPOPLE(NN)
                  ELSE
                     NN = N1+KLGT(KLN)
                     VAL = GHONDO(NN)
                  END IF
                  IF(ABS(VAL).LT.CUTOFF) GO TO 200
                  I1 = LOCI+I
                  I2 = LOCJ+J
                  I3 = LOCK+K
                  I4 = LOCL+L
                  IF (I1 .GE. I2) GO TO 100
                  N = I1
                  I1 = I2
                  I2 = N
  100             IF (I3 .GE. I4) GO TO 110
                  N = I3
                  I3 = I4
                  I4 = N
  110             IF (MIXED) GOTO 180
                  IF (I1-I3) 140,160,180
  140             N = I1
                  I1 = I3
                  I3 = N
                  N = I2
                  I2 = I4
                  I4 = N
                  GO TO 180
  160             IF (I2 .LT. I4) GO TO 140
  180             CONTINUE
C
                  IF (OUT) CALL INTOUT(I1,I2,I3,I4,QQ4,NN,VAL)
                  IF(.NOT.MIXED) THEN
                     IF (I1 .EQ. I2) VAL = VAL*HALF
                     IF (I3 .EQ. I4) VAL = VAL*HALF
                     IF (I1 .EQ. I3 .AND. I2 .EQ. I4) VAL = VAL*HALF
                  END IF
C
               NPACK = ICOUNT
               IPACK = I1
               JPACK = I2
               KPACK = I3
               LPACK = I4
               IF (LABSIZ .EQ. 2) THEN
*I32             LABEL1 = ISHFT( IPACK, 16 ) + JPACK
*I32             LABEL2 = ISHFT( KPACK, 16 ) + LPACK
*I32             IX( 2*NPACK-1 ) = LABEL1
*I32             IX( 2*NPACK   ) = LABEL2
*I64             LABEL = ISHFT( IPACK, 48 ) + ISHFT( JPACK, 32 ) +
*I64 *                   ISHFT( KPACK, 16 ) + LPACK
*I64             IX(NPACK) = LABEL
               ELSE IF (LABSIZ .EQ. 1) THEN
*I32             LABEL = ISHFT( IPACK, 24 ) + ISHFT( JPACK, 16 ) +
*I32 *                   ISHFT( KPACK,  8 ) + LPACK
*I32             IX(NPACK) = LABEL
*I64             IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64               LABEL = ISHFT( IPACK, 24 ) + ISHFT( JPACK, 16 ) +
*I64 *                     ISHFT( KPACK,  8 ) + LPACK
*I64               IX( NPACK/2 ) = IX( NPACK/2 ) + LABEL
*I64             ELSE
*I64               LABEL = ISHFT( IPACK, 56 ) + ISHFT( JPACK, 48 ) +
*I64 *                     ISHFT( KPACK, 40 ) + ISHFT( LPACK, 32 )
*I64               IX( (NPACK/2)+1 ) = LABEL
*I64             END IF
               END IF
C
                  BUFP(ICOUNT) = VAL
                  ICOUNT = ICOUNT+1
                  IF (ICOUNT .GT. NINTMX) THEN
                     NXX = NINTMX
                     CALL PWRIT(IS,BUFP,IX,NXX,NINTMX)
                     ICOUNT = 1
                     NREC = NREC+1
                  END IF
  200          CONTINUE
  220       CONTINUE
  240    CONTINUE
  260 CONTINUE
      RETURN
      END
C
C*MODULE NEOINT  *DECK NEO_SHELLS
      SUBROUTINE NEO_SHELLS(NELEC,ISH,JSH,KSH,LSH,FIRST)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL IANDJ,KANDL,SAME,FIRST
C
      DIMENSION IX(84),IY(84),IZ(84),
     *          JX(84),JY(84),JZ(84),
     *          KX(84),KY(84),KZ(84),
     *          LX(84),LY(84),LZ(84)
C
      PARAMETER (MXSH=5000, MXGSH=30, MXGTOT=20000, MXATM=2000)
C
      COMMON /GOUT  / GPOPLE(768),NORG
      COMMON /INTDEX/ IJGT(784),IJX(784),IJY(784),IJZ(784),IK(784),
     *                KLGT(784),KLX(784),KLY(784),KLZ(784)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /MISC  / IANDJ,KANDL,SAME
      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 /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     +                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     +                NIJ,IJ,KL,IJKL
      COMMON /SHLINF/ GA(MXGSH),CSA(MXGSH),CPA(MXGSH),CDA(MXGSH),
     *                CFA(MXGSH),CGA(MXGSH),CHA(MXGSH),CIA(MXGSH),
     *                GB(MXGSH),CSB(MXGSH),CPB(MXGSH),CDB(MXGSH),
     *                CFB(MXGSH),CGB(MXGSH),CHB(MXGSH),CIB(MXGSH),
     *                GC(MXGSH),CSC(MXGSH),CPC(MXGSH),CDC(MXGSH),
     *                CFC(MXGSH),CGC(MXGSH),CHC(MXGSH),CIC(MXGSH),
     *                GD(MXGSH),CSD(MXGSH),CPD(MXGSH),CDD(MXGSH),
     *                CFD(MXGSH),CGD(MXGSH),CHD(MXGSH),CID(MXGSH),
     *                AX,AY,AZ,BX,BY,BZ,RAB,CX,CY,CZ,DX,DY,DZ,RCD,
     *                NGA,NGB,NGC,NGD
C
      SAVE IGT,JGT,KGT,LGT
      DATA LX /   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,
     *            5,   0,   0,   4,   4,   1,   0,   1,   0,   3,
     *            3,   2,   0,   2,   0,   3,   1,   1,   2,   2,
     *            1,
     *            6,   0,   0,   5,   5,   1,   0,   1,   0,   4,
     *            4,   2,   0,   2,   0,   4,   1,   1,   3,   3,
     *            0,   3,   3,   2,   1,   2,   1,   2/
      DATA KX /   0,   7,   0,   0,  14,   0,   0,   7,   7,   0,
     *           21,   0,   0,  14,  14,   7,   0,   7,   0,   7,
     *           28,   0,   0,  21,  21,   7,   0,   7,   0,  14,
     *           14,   0,  14,   7,   7,
     *           35,   0,   0,  28,  28,   7,   0,   7,   0,  21,
     *           21,  14,   0,  14,   0,  21,   7,   7,  14,  14,
     *            7,
     *           42,   0,   0,  35,  35,   7,   0,   7,   0,  28,
     *           28,  14,   0,  14,   0,  28,   7,   7,  21,  21,
     *            0,  21,  21,  14,   7,  14,   7,  14/
      DATA JX /   0,  49,   0,   0,  98,   0,   0,  49,  49,   0,
     *          147,   0,   0,  98,  98,  49,   0,  49,   0,  49,
     *          196,   0,   0, 147, 147,  49,   0,  49,   0,  98,
     *           98,   0,  98,  49,  49,
     *          245,   0,   0, 196, 196,  49,   0,  49,   0, 147,
     *          147,  98,   0,  98,   0, 147,  49,  49,  98,  98,
     *           49,
     *          294,   0,   0, 245, 245,  49,   0,  49,   0, 196,
     *          196,  98,   0,  98,   0, 196,  49,  49, 147, 147,
     *            0, 147, 147,  98,  49,  98,  49,  98/
      DATA IX /   1, 344,   1,   1, 687,   1,   1, 344, 344,   1,
     *         1030,   1,   1, 687, 687, 344,   1, 344,   1, 344,
     *         1373,   1,   1,1030,1030, 344,   1, 344,   1, 687,
     *          687,   1, 687, 344, 344,
     *         1716,   1,   1,1373,1373, 344,   1, 344,   1,1030,
     *         1030, 687,   1, 687,   1,1030, 344, 344, 687, 687,
     *          344,
     *         2059,   1,   1,1716,1716, 344,   1, 344,   1,1373,
     *         1373, 687,   1, 687,   1,1373, 344, 344,1030,1030,
     *            1,1030,1030, 687, 344, 687, 344, 687/
      DATA LY /   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,
     *            0,   5,   0,   1,   0,   4,   4,   0,   1,   2,
     *            0,   3,   3,   0,   2,   1,   3,   1,   2,   1,
     *            2,
     *            0,   6,   0,   1,   0,   5,   5,   0,   1,   2,
     *            0,   4,   4,   0,   2,   1,   4,   1,   3,   0,
     *            3,   2,   1,   3,   3,   1,   2,   2/
      DATA KY /   0,   0,   7,   0,   0,  14,   0,   7,   0,   7,
     *            0,  21,   0,   7,   0,  14,  14,   0,   7,   7,
     *            0,  28,   0,   7,   0,  21,  21,   0,   7,  14,
     *            0,  14,   7,  14,   7,
     *            0,  35,   0,   7,   0,  28,  28,   0,   7,  14,
     *            0,  21,  21,   0,  14,   7,  21,   7,  14,   7,
     *           14,
     *            0,  42,   0,   7,   0,  35,  35,   0,   7,  14,
     *            0,  28,  28,   0,  14,   7,  28,   7,  21,   0,
     *           21,  14,   7,  21,  21,   7,  14,  14/
      DATA JY /   0,   0,  49,   0,   0,  98,   0,  49,   0,  49,
     *            0, 147,   0,  49,   0,  98,  98,   0,  49,  49,
     *            0, 196,   0,  49,   0, 147, 147,   0,  49,  98,
     *            0,  98,  49,  98,  49,
     *            0, 245,   0,  49,   0, 196, 196,   0,  49,  98,
     *            0, 147, 147,   0,  98,  49, 147,  49,  98,  49,
     *           98,
     *            0, 294,   0,  49,   0, 245, 245,   0,  49,  98,
     *            0, 196, 196,   0,  98,  49, 196,  49, 147,   0,
     *          147,  98,  49, 147, 147,  49,  98,  98/
      DATA IY /   1,   1, 344,   1,   1, 687,   1, 344,   1, 344,
     *            1,1030,   1, 344,   1, 687, 687,   1, 344, 344,
     *            1,1373,   1, 344,   1,1030,1030,   1, 344, 687,
     *            1, 687, 344, 687, 344,
     *            1,1716,   1, 344,   1,1373,1373,   1, 344, 687,
     *            1,1030,1030,   1, 687, 344,1030, 344, 687, 344,
     *          687,
     *            1,2059,   1, 344,   1,1716,1716,   1, 344, 687,
     *            1,1373,1373,   1, 687, 344,1373, 344,1030,   1,
     *         1030, 687, 344,1030,1030, 344, 687, 687/
      DATA LZ /   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,
     *            0,   0,   5,   0,   1,   0,   1,   4,   4,   0,
     *            2,   0,   2,   3,   3,   1,   1,   3,   1,   2,
     *            2,
     *            0,   0,   6,   0,   1,   0,   1,   5,   5,   0,
     *            2,   0,   2,   4,   4,   1,   1,   4,   0,   3,
     *            3,   1,   2,   1,   2,   3,   3,   2/
      DATA KZ /   0,   0,   0,   7,   0,   0,  14,   0,   7,   7,
     *            0,   0,  21,   0,   7,   0,   7,  14,  14,   7,
     *            0,   0,  28,   0,   7,   0,   7,  21,  21,   0,
     *           14,  14,   7,   7,  14,
     *            0,   0,  35,   0,   7,   0,   7,  28,  28,   0,
     *           14,   0,  14,  21,  21,   7,   7,  21,   7,  14,
     *           14,
     *            0,   0,  42,   0,   7,   0,   7,  35,  35,   0,
     *           14,   0,  14,  28,  28,   7,   7,  28,   0,  21,
     *           21,   7,  14,   7,  14,  21,  21,  14/
      DATA JZ /   0,   0,   0,  49,   0,   0,  98,   0,  49,  49,
     *            0,   0, 147,   0,  49,   0,  49,  98,  98,  49,
     *            0,   0, 196,   0,  49,   0,  49, 147, 147,   0,
     *           98,  98,  49,  49,  98,
     *            0,   0, 245,   0,  49,   0,  49, 196, 196,   0,
     *           98,   0,  98, 147, 147,  49,  49, 147,  49,  98,
     *           98,
     *            0,   0, 294,   0,  49,   0,  49, 245, 245,   0,
     *           98,   0,  98, 196, 196,  49,  49, 196,   0, 147,
     *          147,  49,  98,  49,  98, 147, 147,  98/
      DATA IZ /   1,   1,   1, 344,   1,   1, 687,   1, 344, 344,
     *            1,   1,1030,   1, 344,   1, 344, 687, 687, 344,
     *            1,   1,1373,   1, 344,   1, 344,1030,1030,   1,
     *          687, 687, 344, 344, 687,
     *            1,   1,1716,   1, 344,   1, 344,1373,1373,   1,
     *          687,   1, 687,1030,1030, 344, 344,1030, 344, 687,
     *          687,
     *            1,   1,2059,   1, 344,   1, 344,1716,1716,   1,
     *          687,   1, 687,1373,1373, 344, 344,1373,   1,1030,
     *         1030, 344, 687, 344, 687,1030,1030, 687/
C
CECI      DATA FIRST/.TRUE./
C
      IF(FIRST) THEN
         FIRST=.FALSE.
         CALL BASCHK(LMAX)
                       NANGM =  4
         IF(LMAX.EQ.2) NANGM =  6
         IF(LMAX.EQ.3) NANGM = 10
         IF(LMAX.EQ.4) NANGM = 15
         LGT = 1
         KGT = LGT * NANGM
         JGT = KGT * NANGM
         IGT = JGT * NANGM
      END IF
C
C     PREPARE SHELL INFORMATION/FOR HONDO INTEGRATION
C
      IF(NELEC.EQ.2) GO TO 200
C
C     ----- PERMUTE ISH AND JSH SHELLS, FOR THEIR TYPE
C
      IANDJ = ISH .EQ. JSH
      IF (KTYPE(ISH) .LT. KTYPE(JSH)) THEN
         INU = JSH
         JNU = ISH
         NGTI = JGT
         NGTJ = IGT
      ELSE
         INU = ISH
         JNU = JSH
         NGTI = IGT
         NGTJ = JGT
      END IF
C
C     ----- ISHELL
C
      I = KATOM(INU)
      AX = C(1,I)
      AY = C(2,I)
      AZ = C(3,I)
      I1 = KSTART(INU)
      I2 = I1+KNG(INU)-1
      LIT = KTYPE(INU)
      MINI = KMIN(INU)
      MAXI = KMAX(INU)
      LOCI = KLOC(INU)-MINI
      NGA = 0
      DO 140 I = I1,I2
         NGA = NGA+1
         GA(NGA) = EX(I)
         CSA(NGA) = CS(I)
         CPA(NGA) = CP(I)
         CDA(NGA) = CD(I)
         CFA(NGA) = CF(I)
         CGA(NGA) = CG(I)
  140 CONTINUE
C
C     ----- JSHELL
C
      J = KATOM(JNU)
      BX = C(1,J)
      BY = C(2,J)
      BZ = C(3,J)
      J1 = KSTART(JNU)
      J2 = J1+KNG(JNU)-1
      LJT = KTYPE(JNU)
      MINJ = KMIN(JNU)
      MAXJ = KMAX(JNU)
      LOCJ = KLOC(JNU)-MINJ
      NGB = 0
      DO 160 J = J1,J2
         NGB = NGB+1
         GB(NGB) = EX(J)
         CSB(NGB) = CS(J)
         CPB(NGB) = CP(J)
         CDB(NGB) = CD(J)
         CFB(NGB) = CF(J)
         CGB(NGB) = CG(J)
  160 CONTINUE
      RAB = ((AX-BX)*(AX-BX) + (AY-BY)*(AY-BY) + (AZ-BZ)*(AZ-BZ))
C
C     ----- PREPARE INDICES FOR PAIRS OF (I,J) FUNCTIONS
C
      IJ = 0
      JMAX = MAXJ
      DO 190 I = MINI,MAXI
         NX = IX(I)
         NY = IY(I)
         NZ = IZ(I)
         IF (IANDJ) JMAX = I
         DO 180 J = MINJ,JMAX
            IJ = IJ+1
            IJX(IJ) = NX+JX(J)
            IJY(IJ) = NY+JY(J)
            IJZ(IJ) = NZ+JZ(J)
            IJGT(IJ) = NGTI*(I-MINI)+NGTJ*(J-MINJ)+1+NORG
  180    CONTINUE
  190 CONTINUE
      RETURN
C     ******
C
C        K AND L SHELL
C
  200 CONTINUE
      KANDL = KSH .EQ. LSH
      SAME = ISH .EQ. KSH .AND. JSH .EQ. LSH
C
C     ----- PERMUTE KSH AND LSH SHELLS, FOR THEIR TYPE
C
      IF (KTYPE(KSH) .LT. KTYPE(LSH)) THEN
         KNU = LSH
         LNU = KSH
         NGTK = LGT
         NGTL = KGT
      ELSE
         KNU = KSH
         LNU = LSH
         NGTK = KGT
         NGTL = LGT
      END IF
C
C     ----- K SHELL
C
      K = KATOM(KNU)
      CX = C(1,K)
      CY = C(2,K)
      CZ = C(3,K)
      K1 = KSTART(KNU)
      K2 = K1+KNG(KNU)-1
      LKT = KTYPE(KNU)
      MINK = KMIN(KNU)
      MAXK = KMAX(KNU)
      LOCK = KLOC(KNU)-MINK
      NGC = 0
      DO 260 K = K1,K2
         NGC = NGC+1
         GC(NGC) = EX(K)
         CSC(NGC) = CS(K)
         CPC(NGC) = CP(K)
         CDC(NGC) = CD(K)
         CFC(NGC) = CF(K)
         CGC(NGC) = CG(K)
  260 CONTINUE
C
C     ----- LSHELL
C
      L = KATOM(LNU)
      DX = C(1,L)
      DY = C(2,L)
      DZ = C(3,L)
      L1 = KSTART(LNU)
      L2 = L1+KNG(LNU)-1
      LLT = KTYPE(LNU)
      MINL = KMIN(LNU)
      MAXL = KMAX(LNU)
      LOCL = KLOC(LNU)-MINL
      NGD = 0
      DO 280 L = L1,L2
         NGD = NGD+1
         GD(NGD) = EX(L)
         CSD(NGD) = CS(L)
         CPD(NGD) = CP(L)
         CDD(NGD) = CD(L)
         CFD(NGD) = CF(L)
         CGD(NGD) = CG(L)
  280 CONTINUE
      NROOTS = (LIT+LJT+LKT+LLT-2)/2
      RCD = ((CX-DX)*(CX-DX) + (CY-DY)*(CY-DY) + (CZ-DZ)*(CZ-DZ))
C
C     ----- PREPARE INDICES FOR PAIRS OF (K,L) FUNCTIONS
C
      KL = 0
      LMAX = MAXL
      DO 310 K = MINK,MAXK
         NX = KX(K)
         NY = KY(K)
         NZ = KZ(K)
         IF (KANDL) LMAX = K
         DO 300 L = MINL,LMAX
            KL = KL+1
            KLX(KL) = NX+LX(L)
            KLY(KL) = NY+LY(L)
            KLZ(KL) = NZ+LZ(L)
            KLGT(KL) = NGTK*(K-MINK)+NGTL*(L-MINL)
  300    CONTINUE
  310 CONTINUE
      MAX = KL
      DO 320 I = 1,IJ
      IF (SAME) MAX = I
  320 IK(I) = MAX
      IJKL = IJ*KL
      IF (SAME) IJKL = IJ*(IJ+1)/2
      RETURN
      END
C
C*MODULE NEOINT  *DECK NEO_ZQOUT
      SUBROUTINE NEO_ZQOUT(GHONDO,MIXED)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION GHONDO(*)
C
      LOGICAL IANDJ,KANDL,SAME,MIXED
C
      COMMON /INTDEX/ IJGT(784),IJX(784),IJY(784),IJZ(784),IK(784),
     *                KLGT(784),KLX(784),KLY(784),KLZ(784)
      COMMON /MISC  / IANDJ,KANDL,SAME
      COMMON /SHLNOS/ QQ4,LIT,LJT,LKT,LLT,LOCI,LOCJ,LOCK,LOCL,
     *                MINI,MINJ,MINK,MINL,MAXI,MAXJ,MAXK,MAXL,
     *                NIJX,IJ,KL,IJKL
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- ZERO HONDO CONVENTIONAL INTEGRAL OUTPUT REGION -----
C
      IJN = 0
      JMAX = MAXJ
      DO 260 I = MINI,MAXI
         IF (IANDJ) JMAX = I
         DO 240 J = MINJ,JMAX
            IJN = IJN+1
            N1 = IJGT(IJN)
            LMAX = MAXL
            KLN = 0
            DO 220 K =  MINK,MAXK
               IF (KANDL) LMAX = K
               DO 200 L = MINL,LMAX
                  KLN = KLN+1
                  IF(.NOT.MIXED) THEN
                     IF (SAME .AND. KLN .GT. IJN) GO TO 240
                  END IF
                  NN = N1+KLGT(KLN)
                  GHONDO(NN) = ZERO
  200          CONTINUE
  220       CONTINUE
  240    CONTINUE
  260 CONTINUE
      RETURN
      END
