C 08 AUG 07 - JHS - ADD PBFLAB TO COMMON BLOCK NUCMON
C 07 AUG 07 - JHS - REPLACE BFLAB LABEL FOR PROTONS WITH PFLAB
C 07 AUG 07 - JHS - ADDED QUADRUPLE ZETA BASIS SETS
C 25 OCT 06 - SPW - CHANGES FOR NEO VERSION 2
C 25 MAR 06 - SPW - NEW MODULE FOR NEO BASIS SETS
C
C*MODULE NEOBAS  *DECK NBASX
      SUBROUTINE NBASX
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
CSPW  LOGICAL USRDEX
C
      PARAMETER (MXGTOT=20000)
C
CSPW  COMMON /NBSEXP/ EXNB(MAXNUP,MAXQNT),CONNB(MAXNUP,MAXQNT),
CSPW *                ICONNB(MAXNUP,MAXQNT),USRDEX
      COMMON /FMCOM / X(1)
C
      CALL VALFM(LOADFM)
      LS   = LOADFM + 1
      LP   = LS + MXGTOT
      LD   = LP + MXGTOT
      LF   = LD + MXGTOT
      LG   = LF + MXGTOT
      LAST = LG + MXGTOT
      NEED = LAST - LS
      CALL GETFM(NEED)
C
      CALL NUCBAS(X(LS),X(LP),X(LD),X(LF),X(LG),MXGTOT)
      CALL RETFM(NEED)
C
      RETURN
      END
C
C*MODULE NEOBAS  *DECK NUCBAS
      SUBROUTINE NUCBAS(CSINP,CPINP,CDINP,CFINP,CGINP,MAXGAU)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
CSPW  LOGICAL USRDEX
C
      PARAMETER (MXSH=5000, MXGSH=30, MXGTOT=20000, MXAO=8192,
     *          MXATM=2000, MXRT=100)
C
      CHARACTER*4 ATMSYM(MXATM)
      CHARACTER*8 QNUN,QNN
      CHARACTER*10 PBFLAB
C
      DOUBLE PRECISION LINE,LINEV,LABEL(28)
C
      DIMENSION CSINP(MAXGAU),CPINP(MAXGAU),CDINP(MAXGAU),
     *          CFINP(MAXGAU),CGINP(MAXGAU)
      DIMENSION INTYP(MXSH),NS(MXATM),KS(MXATM),
     *          NBFS(27),MINF(27),MAXF(27),NANGM(27)
C
CSPW  COMMON /NBSEXP/ EXNB(MAXNUP,MAXQNT),CONNB(MAXNUP,MAXQNT),
CSPW *                ICONNB(MAXNUP,MAXQNT),USRDEX
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      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 /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 /NEOSHL/ NGAUSS,NGAUSE,NGAUSN,NTSHEL,NNUCSH
      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      COMMON /NUCMOM/ BASNUC
      COMMON /NUCMON/ QNUN(20),QNN(20),PBFLAB(MXAO)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),A(MXATM),B(MXATM),BFLAB(MXAO)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TRANSF/ XSMAL,YSMAL,ZSMAL,XNEW,YNEW,ZNEW,XP,YP,ZP
C
      PARAMETER (ZERO=0.0D+00, TM6=1.0D-06, TM10=1.0D-10)
C
      DATA SZSNB,DZSNB,TZSNB,QZSNB/8HSZSNB   ,8HDZSNB   ,8HTZSNB   ,
     *     8HQZSNB   /
      DATA SZSPNB,DZSPNB,TZSPNB,QZSPNB/8HSZSPNB  ,8HDZSPNB  ,8HTZSPNB  ,
     *     8HQZSPNB  /
      DATA SZSPDN,DZSPDN,TZSPDN,QZSPDD/8HSZSPDN  ,8HDZSPDN  ,
     *     8HTZSPDN  ,8HQZSPDD  /
      DATA GENNBS/8HGENNBS  /
      DATA READIN/8HREADIN  /
C
C           THIS ROUTINE IS PARTIALLY SET TO READ F AND N ORBITALS
C           BUT THESE ARE NOT IMPLEMETED IN THE INTEGRAL
C           SECTIONS OF THIS PROGRAM, AS OF YET.
C
      DATA LABEL/8HS       ,8HP       ,8HD       ,8HF       ,
     *           8HG       ,8HL       ,8HM       ,8HN       ,
     *           8H1S      ,8H2S      ,8H2P      ,8H2SP     ,
     *           8H3S      ,8H3P      ,8H3D      ,8H3SP     ,8H3SPD    ,
     *           8H4S      ,8H4P      ,8H4D      ,8H4SP     ,8H4SPD    ,
     *           8H5S      ,8H5P      ,8H5D      ,8H5SP     ,8H5SPD    ,
     *           8HSV      /
C
      DATA NBFS/ 1, 3, 6,10,15,     4,10,20,
     *           1, 1, 3, 4,        1, 3, 6, 4,10,
     *           1, 3, 6, 4,10,     1, 3, 6, 4,10/
      DATA MINF/ 1, 2, 5,11,21,     1, 1, 1,
     *           1, 1, 2, 1,        1, 2, 5, 1, 1,
     *           1, 2, 5, 1, 1,     1, 2, 5, 1, 1/
      DATA MAXF/ 1, 4,10,20,35,     4,10,20,
     *           1, 1, 4, 4,        1, 4,10, 4,10,
     *           1, 4,10, 4,10,     1, 4,10, 4,10/
      DATA NANGM/1, 2, 3, 4, 5,     2, 3, 4,
     *           1, 1, 2, 2,        1, 2, 3, 2, 3,
     *           1, 2, 3, 2, 3,     1, 2, 3, 2, 3/
      DATA LINE,LINEV /8H * * * *,8H       */
C
      IGAUSS = 0
C
C     READ IN CENTERS AND NUCLEAR BASIS SETS GROUPED IN SHELLS
C
      NNUCSH = 0
      NTSHEL = NSHELL
      LOC = 0
      NGAUSE = NGAUSS
      X = ZERO
      Y = ZERO
      Z = ZERO
      NQN = 0
C
C     --- LOOP OVER QUANTUM NUCLEI BASIS SET CENTRES ---
C
      IF (MASWRK) THEN
         WRITE (IW,7080)
      END IF
C
      NTOTCN = NUNIQN + NAUXNB
C
      DO 5000 II = 1,NTOTCN
C
         NQN = NQN + 1
C
         IF(II.LE.NUNIQN) THEN
            IAT = IUNIQN(II)
         ELSE
            IAT = IAUXNB(II-NUNIQN)
         END IF
         NUCTYP = IUNIQT(II)
         QNN(II) = QNUN(II)
         IQMNUC(NQN) = IAT
         IQNTYP(NQN) = NUCTYP
C
         X = C(1,IAT)
         Y = C(2,IAT)
         Z = C(3,IAT)
         NS(NQN) = 0
         KS(NQN) = NNUCSH+1
C
         IF(IGAUSS .GT. MXGSH) THEN
            IF (MASWRK) WRITE (IW,9010) MXGSH
            CALL ABRT
         END IF
C
C        --- SET UP NUCLEAR AB INITIO BASIS SET -BASNUC- ---
C
         IERR1=0
         IERR2=0
C
C        --- SZSNB (SINGLE-ZETA S) ---
C
         IF (BASNUC .EQ. SZSNB) THEN
C
           NN = 1
           CALL SZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *               INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *               IAT,NTSHEL,NNUCSH,NQN)
C
C        --- DZSNB (DOUBLE-ZETA S) ---
C
         ELSE IF (BASNUC .EQ. DZSNB) THEN
C
           NN = 1
           CALL DZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *               INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *               IAT,NTSHEL,NNUCSH,NQN)
C
C        --- TZSNB (TRIPLE-ZETA S) ---
C
         ELSE IF (BASNUC .EQ. TZSNB) THEN
C
            NN=1
            CALL TZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *                IAT,NTSHEL,NNUCSH,NQN)
C
C        --- QZSNB (QUADRUPLE-ZETA S) ---
C
         ELSE IF (BASNUC .EQ. QZSNB) THEN
C
           NN = 1
           CALL QZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *               INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *               IAT,NTSHEL,NNUCSH,NQN)
C
C        --- SZSPNB (SINGLE-ZETA SP) ---
C
         ELSE IF (BASNUC .EQ. SZSPNB) THEN
C
           NN = 2
           CALL SZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *               INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *               IAT,NTSHEL,NNUCSH,NQN)
C
C        --- DZSPNB (DOUBLE-ZETA SP) ---
C
         ELSE IF (BASNUC .EQ. DZSPNB) THEN
C
            NN = 2
            CALL DZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *                IAT,NTSHEL,NNUCSH,NQN)
C
C        --- TZSPNB (TRIPLE-ZETA SP) ---
C
         ELSE IF(BASNUC .EQ. TZSPNB) THEN
C
            NN = 2
            CALL TZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *                IAT,NTSHEL,NNUCSH,NQN)
C
C        --- QZSPNB (QUADRUPLE-ZETA SP) ---
C
         ELSE IF (BASNUC .EQ. QZSPNB) THEN
C
            NN = 2
            CALL QZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *                IAT,NTSHEL,NNUCSH,NQN)
C
C        --- SZSPDN (SINGLE-ZETA SPD) ---
C
         ELSE IF (BASNUC .EQ. SZSPDN) THEN
C
           NN = 3
           CALL SZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *               INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *               IAT,NTSHEL,NNUCSH,NQN)
C
C        --- DZSPDN (DOUBLE-ZETA SPD) ---
C
         ELSE IF (BASNUC .EQ. DZSPDN) THEN
C
            NN = 3
            CALL DZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *                IAT,NTSHEL,NNUCSH,NQN)
C
C        --- TZSPDN (TRIPLE-ZETA SPD) ---
C
         ELSE IF(BASNUC .EQ. TZSPDN) THEN
C
            NN = 3
            CALL TZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *                IAT,NTSHEL,NNUCSH,NQN)
C
C        --- QZSPDD (QUADRUPLE-ZETA SP  DOUBLE-ZETA D) ---
C
         ELSE IF (BASNUC .EQ. QZSPDD) THEN
C
            NN = 3
            CALL QZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *                IAT,NTSHEL,NNUCSH,NQN)
C
C        --- COMPLETELY USER DEFINED BASIS SET ---
C
         ELSE IF(BASNUC .EQ. GENNBS) THEN
C
            WRITE(IW,*)'USER DEFINED NUCLEAR BASIS OPTION NOT ',
     *                 'IMPLEMENTED YET'
            CALL ABRT
C
C        --- READ IN GENRAL BASIS SET FROM EXTERNAL FILE ---
C
         ELSE IF (BASNUC .EQ. READIN) THEN
C
            CALL READNB(IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                  INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,
     *                  IAT,NTSHEL,NNUCSH,NQN)
C
         ELSE
            CALL ABRT
         END IF
C
C        --- CHECK FOR ERRORS IN BASIS ---
C
         IF(IERR1.NE.0) GO TO 1560
         IF(IERR2.NE.0) GO TO 1600
C
C 800    CONTINUE
C
C        --- GENERATE EQUIVALENT QUANTUM NUCLEAR CENTERS ---
C
         JJA = IAT
         XS = X
         YS = Y
         ZS = Z
         CALL LOCAL(X,Y,Z,XS,YS,ZS)
         XSMAL = XS
         YSMAL = YS
         ZSMAL = ZS
         NQN0 = NQN
C
         DO 820 IT = 1,NT
            IF(IT .EQ. 1) GO TO 820
            NN = 9*(IT-1)
            CALL TRANS(NN)
            CALL ROT
C
            DO IDX = 1,NQN
               IIA = IQMNUC(IDX)
               TEST = (XP-C(1,IIA))**2+(YP-C(2,IIA))**2+(ZP-C(3,IIA))**2
               IF(TEST .LE. TM6) GO TO 820
            END DO
C
            DO IIA = JJA+1,NAT
               TEST = (XP-C(1,IIA))**2+(YP-C(2,IIA))**2+(ZP-C(3,IIA))**2
               IF(TEST .LE. TM10) GO TO 810
            END DO
  810       CONTINUE
            JJA = IIA
C
            NQN = NQN + 1
            IF(NQN.GT.20) GO TO 1540
C
            IQMNUC(NQN) = IIA
            IQNTYP(NQN) = NUCTYP
            QNN(NQN) = QNUN(II)
C
            NS(NQN) = NS(NQN0)
            KS(NQN) = KS(NQN-1) + NS(NQN0)
            NNUCSH = NNUCSH + NS(NQN)
            NTSHEL = NTSHEL + NS(NQN)
            NS1 = NS(NQN)
            J = NSHELL + KS(NQN0)
            JJ = NSHELL + KS(NQN)
C
            CALL ATOMSS(NBFS,IIA,LOC,NS1,KMIN(J),KMAX(J),KSTART(J),
     *                  KTYPE(J),KNG(J),INTYP(J),KATOM(JJ),
     *                  KMIN(JJ),KMAX(JJ),KSTART(JJ),KTYPE(JJ),KNG(JJ),
     *                  INTYP(JJ),KLOC(JJ))
  820    CONTINUE
C
 5000 CONTINUE
C
C     --- ALL QUANTUM NUCLEI LOOPED OVER ---
C
      NQMNUC = NQN - NTAUXB
C
C     --- CHECK BASIS FUNCTIONS DEFINED ---
C
      IF (NGAUSS.EQ.0  .OR.  NNUCSH.EQ.0) THEN
         IF (MASWRK) THEN
            WRITE(IW,*) 'ERROR!  THERE ARE NO BASIS FUNCTIONS DEFINED.'
            WRITE(IW,*) 'CHECK YOUR $BASIS/$DATA INPUT FOR ERRORS.'
         END IF
         CALL ABRT
         STOP
      END IF
C
C     ----- FORM TRANSFORMATION TABLES FOR NUCLEAR SHELLS -----
C
      DO 1050 IQN = 1,NQN
C
         IAT = IQMNUC(IQN)
         NS1 = NSHELL + KS(IQN)-1
         NS2 = NS(IQN)
         X = C(1,IAT)
         Y = C(2,IAT)
         Z = C(3,IAT)
         XS = X
         YS = Y
         ZS = Z
         CALL LOCAL(X,Y,Z,XS,YS,ZS)
         XSMAL = XS
         YSMAL = YS
         ZSMAL = ZS
C
         DO 1040 IT = 1,NT
            NN = 9*(IT-1)
            CALL TRANS(NN)
            CALL ROT
            ICTR = -2**20
            DO 1010 I = 1,NQN
               JAT = IQMNUC(I)
               TEST = (XP-C(1,JAT))**2+(YP-C(2,JAT))**2+(ZP-C(3,JAT))**2
               IF(TEST .GT. TM10) GO TO 1010
               ICTR = I
               GO TO 1020
 1010       CONTINUE
 1020       NS3 = KS(ICTR) - 1
            DO 1030 ISH = 1,NS2
               MAPSHL(NS1+ISH,IT) = NSHELL + NS3 + ISH
 1030       CONTINUE
 1040    CONTINUE
C
 1050 CONTINUE
C
C     ----- PRINT NUCLEAR BASIS SET CENTERS AND SHELL DESCRIPTIONS -----
C
      IF (MASWRK) THEN
         WRITE (IW,9050)
         DO 1100 IQN = 1,NQN
            IAT = IQMNUC(IQN)
            WRITE (IW,9060) QNN(IQN),B(IAT),ZAN(IAT),
     *                      C(1,IAT),C(2,IAT),C(3,IAT)
 1100    CONTINUE
      END IF
C
C     ----- PRINT OUT THE BASIS FUNCTIONS USED -----
C
      IF (MASWRK) THEN
         WRITE (IW,9080)
         IF(NORMP .NE. 1) WRITE (IW,9090)
         IF(NORMF .NE. 1) WRITE (IW,9100)
         WRITE (IW,9110)
         MLDNDA=0
         DO 1420 IQN = 1,NQN
            IAT = IQMNUC(IQN)
            DO 1260 IT = 1,NT
               IF(MAPCTR(IAT,IT) .GT. IAT) GO TO 1420
 1260       CONTINUE
            WRITE (IW,9120) QNN(IQN),B(IAT)
            NS1 = NSHELL + KS(IQN)
            NS2 = NS1+NS(IQN)-1
            MLDNDA=MLDNDA+1
            DO 1270 ISH=1,5
 1270       CONTINUE
            DO 1400 ISH = NS1,NS2
               ISHN = ISH - NSHELL
               WRITE (IW,9130)
               I1 = KSTART(ISH)
               I2 = I1+KNG(ISH)-1
               ITYP = INTYP(ISH)
               DO 1380 IG = I1,I2
                  IGN = IG - KSTART(NSHELL)
C
                  GO TO (1280,1290,1300,1310,
     *                   1315,1320,1330,1340,
     *                   1280,     1280,1290,1320,
     *                   1280,1290,1300,1320,1330,
     *                   1280,1290,1300,1320,1330,
     *                   1280,1290,1300,1320,1330), ITYP
C
 1280             CONTINUE
                  C2=CSINP(IG)
                  WRITE(IW,9140) ISHN,LABEL(ITYP),IGN,EX(IG),C2
                  GO TO 1380
C
 1290             CONTINUE
                  C2=CPINP(IG)
                  WRITE(IW,9140) ISHN,LABEL(ITYP),IGN,EX(IG),C2
                  GO TO 1380
C
 1300             CONTINUE
                  C2=CDINP(IG)
                  WRITE(IW,9140) ISHN,LABEL(ITYP),IGN,EX(IG),C2
                  GO TO 1380
C
 1310             CONTINUE
                  C2=CFINP(IG)
                  WRITE(IW,9140) ISHN,LABEL(ITYP),IGN,EX(IG),C2
                  GO TO 1380
C
 1315             CONTINUE
                  C2=CGINP(IG)
                  WRITE(IW,9140) ISHN,LABEL(ITYP),IGN,EX(IG),C2
                  GO TO 1380
C
 1320             CONTINUE
                  C2=CSINP(IG)
                  C4=CPINP(IG)
                  WRITE (IW,9140) ISHN,LABEL(ITYP),IGN,EX(IG),C2,C4
                  GO TO 1380
C
 1330             CONTINUE
                  C2=CSINP(IG)
                  C4=CPINP(IG)
                  C6=CDINP(IG)
                  WRITE (IW,9140) ISHN,LABEL(ITYP),IGN,EX(IG),C2,C4,C6
                  GO TO 1380
C
 1340             CONTINUE
                  C2=CSINP(IG)
                  C4=CPINP(IG)
                  C6=CDINP(IG)
                  C8=CFINP(IG)
                  WRITE(IW,9140) ISHN,LABEL(ITYP),IGN,EX(IG),C2,C4,C6,C8
C
 1380          CONTINUE
 1400       CONTINUE
 1420    CONTINUE
C
C        ----- OPTIONAL PRINTOUT OF BASIS SYMMETRY INFO -----
C
         IF(NPRINT .NE. 1) GO TO 1530
         WRITE (IW,9150)
         WRITE (IW,9160)
         IMAX = 0
 1440    IMIN = IMAX+1
         IMAX = IMAX+15
         IF(IMAX .GT. NT) IMAX = NT
         IMAX1 = IMAX+1
         WRITE (IW,9170)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9190) LINEV,(I,I = IMIN,IMAX)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINE ,I = IMIN,IMAX1)
         DO 1460 IAT = 1,NAT
            WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
            WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
            WRITE (IW,9200) IAT,(MAPCTR(IAT,I),I = IMIN,IMAX)
            WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
            WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
            WRITE (IW,9180) (LINE ,I = IMIN,IMAX1)
 1460    CONTINUE
         IF(IMAX .LT. NT) GO TO 1440
         WRITE (IW,9150)
         WRITE (IW,9210)
         IMAX = 0
 1480    IMIN = IMAX+1
         IMAX = IMAX+15
         IF(IMAX .GT. NT) IMAX = NT
         IMAX1 = IMAX+1
         WRITE (IW,9170)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9190) LINEV,(I,I = IMIN,IMAX)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINE ,I = IMIN,IMAX1)
         DO 1500 ISH = 1,NSHELL
            WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
            WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
            WRITE (IW,9200) ISH,(MAPSHL(ISH,I),I = IMIN,IMAX)
            WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
            WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
            WRITE (IW,9180) (LINE ,I = IMIN,IMAX1)
 1500    CONTINUE
         IF(IMAX .LT. NT) GO TO 1480
         WRITE (IW,9150)
         WRITE (IW,9220)
         IMAX = 0
 1520    IMIN = IMAX+1
         IMAX = IMAX+15
         IF(IMAX .GT. NT) IMAX = NT
         IMAX1 = IMAX+1
         WRITE (IW,9170)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9190) LINEV,(I,I = IMIN,IMAX)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINE ,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9190) LINEV,(INVT(I),I = IMIN,IMAX)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINEV,I = IMIN,IMAX1)
         WRITE (IW,9180) (LINE ,I = IMIN,IMAX1)
         IF(IMAX .LT. NT) GO TO 1520
      END IF
 1530 CONTINUE
C
C     ---- KILL JOB IF DIMENSIONS HAVE BEEN EXCEEDED ----
C
 1540 CONTINUE
      IF(NAT .LE. MXATM) GO TO 1560
      IF (MASWRK) WRITE (IW,9230) MXATM
      CALL ABRT
C
 1560 IF(NTSHEL .LE. MXSH) GO TO 1580
      IF (MASWRK) WRITE (IW,9240) MXSH
      CALL ABRT
C
 1580 IF(LOC .LE. MXAO) GO TO 1600
      IF (MASWRK) WRITE (IW,9250) MXAO
      CALL ABRT
C
 1600 IF(NGAUSS .LE. MXGTOT) GO TO 1620
      IF (MASWRK) WRITE (IW,9260) MXGTOT
      CALL ABRT
C
C     ESTABLISH NUMBER OF QUANTUM NUCLEI, MULTIPLICITY, ETC.
C     LATER, NQMT MAY BE RESET TO VALUES SMALLER THAN NUM.
C
 1620 CONTINUE
      NUMNB  = LOC
      NGAUSN = NGAUSS - NGAUSE
      NNA = (NQMNUC + NUMULT -1)/2
      NNB = (NQMNUC - NUMULT +1)/2
      IF (MASWRK) THEN
         WRITE(IW,9270) NNUCSH,NUMNB
         WRITE(IW,9272) NUMULT,NNA,NNB,NQMNUC
      END IF
C
C     ----- CORRECT CLASSICAL NUCLEAR REPULSION ENERGY -----
C
      DO IAT = 1,NAT
         DO IQN = 1,NQMNUC
            IF (IAT.EQ.IQMNUC(IQN)) THEN
               ZAN(IAT) = ZERO
            END IF
         END DO
      END DO
      DO IAT = 1,NAT
         DO IAUX = 1,NTAUXB
            IF (IAT.EQ.IAUXNB(IAUX)) THEN
               ZAN(IAT) = ZERO
            END IF
         END DO
      END DO
      ENUCR  = ENUC(NAT,ZAN,C)
      IF (MASWRK) THEN
         WRITE(IW,9280) ENUCR
      END IF
C
      CALL NUBLAB(1,ATMSYM)
C
      RETURN
C
 7080 FORMAT(/5X,'NUCLEAR BASIS SET CENTERS'/5X,25(1H-))
 9010 FORMAT(1X,'*** ERROR, NO MORE THAN',I6,' GAUSSIANS PER SHELL')
 9050 FORMAT(/1X,'NUC ',6X,'ATOMIC',22X,'COORDINATES (BOHR)'/
     *         11X,'CHARGE',9X,'X',19X,'Y',19X,'Z')
 9060 FORMAT(1X,A8,A2,F5.1,F17.10,2F20.10)
 9080 FORMAT(/5X,'NUCLEAR BASIS SET'/5X,17(1H-))
 9090 FORMAT(1X,'THE CONTRACTED PRIMITIVE FUNCTIONS HAVE BEEN',
     *          ' UNNORMALIZED')
 9100 FORMAT(1X,'THE CONTRACTED BASIS FUNCTIONS ARE NOW NORMALIZED',
     *          ' TO UNITY')
 9110 FORMAT(/1X,'SHELL TYPE PRIMITIVE    EXPONENT',
     *       10X,'CONTRACTION COEFFICIENTS')
 9120 FORMAT(/1X,A8,A2)
 9130 FORMAT(1X)
 9140 FORMAT(1X,I3,3X,A2,I4,F20.7,4F20.12)
 9150 FORMAT('1')
 9160 FORMAT(/' **** TRANSFORMATION TABLE OF ATOMS *****'/
     *       30X,15H ROWS ARE ATOMS/
     *       30X,'COLUMNS ARE SYMMETRY OPERATIONS')
 9180 FORMAT(1X,16A8)
 9190 FORMAT(1X,A8,15(2X,I2,3X,1H*))
 9170 FORMAT(/)
 9200 FORMAT(1X,16(2X,I2,3X,1H*))
 9210 FORMAT(/1X,5H*****,31H TRANSFORMATION TABLE OF SHELLS, 5H*****/
     *        30X,16H ROWS ARE SHELLS/
     *        30X,32H COLUMNS ARE SYMMETRY OPERATIONS)
 9220 FORMAT(/1X,'***** INVERSE TRANSFORMATIONS *****')
 9230 FORMAT(1X,'*** ERROR, NO MORE THAN',I6,' TOTAL ATOMS')
 9240 FORMAT(1X,'*** ERROR, NO MORE THAN',I6,' TOTAL SHELLS')
 9250 FORMAT(1X,'*** ERROR, NO MORE THAN',I6,' TOTAL BASIS FUNCTIONS')
 9260 FORMAT(1X,'*** ERROR, NO MORE THAN',I6,' TOTAL GAUSSIANS')
 9270 FORMAT(/1X,'TOTAL NUMBER OF NUCLEAR BASIS SET SHELLS     =',I5/
     *        1X,'NUMBER OF CARTESIAN GAUSSIAN BASIS FUNCTIONS =',I5)
 9272 FORMAT( 1X,'NUCLEAR SPIN MULTIPLICITY                    =',I5/
     *        1X,'NUMBER OF OCCUPIED NUCLEAR ORBITALS (ALPHA)  =',I5/
     *        1X,'NUMBER OF OCCUPIED NUCLEAR ORBITALS (BETA )  =',I5/
     *        1X,'TOTAL NUMBER OF QUANTUM NUCLEI               =',I5)
 9280 FORMAT(1X,'THE CORRECTED CLASSICAL NUCLEAR REPULSION ENERGY IS',
     *       F20.10)
      END
C
C*MODULE NEOBAS  *DECK DZNB
      SUBROUTINE DZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,
     *                NUCTYP,IAT,NTSHEL,NNUCSH,NQN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION CSINP(*),CPINP(*),CDINP(*),
     *          INTYP(*),NANGM(*),NBFS(*),MINF(*),MAXF(*),NS(*)
      DIMENSION EEX(30),CCS(30),CCP(30),CCD(30),
     *          ITYPES(11),NGBFS(11)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      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
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, PT5=0.5D+00, PT75=0.75D+00)
      PARAMETER (PI32=5.56832799683170D+00, TOL=1.0D-10)
C
C     ----- SET UP THE "DOUBLE ZETA" NUCLEAR BASIS SET -----
C
      DO 100 I = 1,30
         EEX(I) = ZERO
         CCS(I) = ZERO
         CCP(I) = ZERO
         CCD(I) = ZERO
  100 CONTINUE
C
C     ----- PROTON -----
C
      IF (NUCTYP.EQ.1) THEN
C
         CALL NDZV01(EEX,CCS,CCP,CCD,NN)
C
C     ----- DEUTERIUM NUCLEUS -----
C
      ELSE IF (NUCTYP.EQ.2) THEN
C
         CALL ABRT
C
C     ----- TRITIUM NUCLEUS -----
C
      ELSE IF (NUCTYP.EQ.3) THEN
C
         CALL ABRT
C
      ELSE
         CALL ABRT
      END IF
C
C     ----- NOW SET UP EACH SHELL -----
C
      CALL NDZSHL(NN,MXPASS,ITYPES,NGBFS)
C
      IGX = 0
      IPASS = 0
  310 CONTINUE
      IPASS = IPASS + 1
      IF(IPASS.GT.MXPASS) RETURN
C                         ******
      IGAUSS = NGBFS(IPASS)
      ITYP = ITYPES(IPASS)
C
      NNUCSH = NNUCSH+1
      NTSHEL = NTSHEL+1
      IF(NTSHEL.GT.MXSH) THEN
         IERR1=1
         RETURN
      END IF
      NS(NQN) = NS(NQN)+1
      KMIN(NTSHEL) = MINF(ITYP)
      KMAX(NTSHEL) = MAXF(ITYP)
      KSTART(NTSHEL) = NGAUSS+1
      KATOM(NTSHEL) = IAT
      KTYPE(NTSHEL) = NANGM(ITYP)
      INTYP(NTSHEL) = ITYP
      KNG(NTSHEL) = IGAUSS
      KLOC(NTSHEL) = LOC+1
C
      NGAUSS = NGAUSS+IGAUSS
      IF(NGAUSS.GT.MXGTOT) THEN
         IERR2=1
         RETURN
      END IF
      LOC = LOC+NBFS(ITYP)
      K1 = KSTART(NTSHEL)
      K2 = K1+KNG(NTSHEL)-1
      DO 720 I = 1,IGAUSS
         K = K1+I-1
         EX(K)    = EEX(IGX+I)
         CSINP(K) = CCS(IGX+I)
         CPINP(K) = CCP(IGX+I)
         CDINP(K) = CCD(IGX+I)
         CS(K) = CSINP(K)
         CP(K) = CPINP(K)
         CD(K) = CDINP(K)
  720 CONTINUE
      IGX = IGX + IGAUSS
C
C     ----- ALWAYS UNNORMALIZE PRIMITIVES -----
C
      DO 740 K = K1,K2
         EE = EX(K)+EX(K)
         FACS = PI32/(EE*SQRT(EE))
         FACP = PT5*FACS/EE
         FACD = PT75*FACS/(EE*EE)
         CS(K) = CS(K)/SQRT(FACS)
         CP(K) = CP(K)/SQRT(FACP)
         CD(K) = CD(K)/SQRT(FACD)
  740 CONTINUE
C
C     ----- IF(NORMF.EQ.0) NORMALIZE BASIS FUNCTIONS. -----
C
      IF (NORMF .EQ. 1) GO TO 820
      FACS = ZERO
      FACP = ZERO
      FACD = ZERO
      DO 780 IG = K1,K2
         DO 770 JG = K1,IG
            EE = EX(IG)+EX(JG)
            FAC = EE*SQRT(EE)
            DUMS = CS(IG)*CS(JG)/FAC
            DUMP = PT5*CP(IG)*CP(JG)/(EE*FAC)
            DUMD = PT75*CD(IG)*CD(JG)/(EE*EE*FAC)
            IF (IG .NE. JG) THEN
               DUMS = DUMS+DUMS
               DUMP = DUMP+DUMP
               DUMD = DUMD+DUMD
            END IF
            FACS = FACS+DUMS
            FACP = FACP+DUMP
            FACD = FACD+DUMD
  770    CONTINUE
  780 CONTINUE
      IF (FACS .GT. TOL) FACS = ONE/SQRT(FACS*PI32)
      IF (FACP .GT. TOL) FACP = ONE/SQRT(FACP*PI32)
      IF (FACD .GT. TOL) FACD = ONE/SQRT(FACD*PI32)
      DO 800 IG = K1,K2
         CS(IG) = CS(IG) * FACS
         CP(IG) = CP(IG) * FACP
         CD(IG) = CD(IG) * FACD
  800 CONTINUE
C
C        AND BRANCH BACK FOR THE NEXT SHELL
C
  820 CONTINUE
      GO TO 310
      END
C
C*MODULE NEOBAS  *DECK NDZSHL
      SUBROUTINE NDZSHL(NN,MXPASS,ITYPES,NGBFS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION ITYPES(11),NGBFS(11)
      DIMENSION KINDS(36),NPASS(8),ITYPS(11,8),NGS(11,8)
C
C     KINDS ARE S, SP, SPD
C     NPASS TELLS HOW MANY SHELLS ARE IN EACH BASIS
C
      DATA KINDS/1,  2,   3,   4,  5,  31*0/
      DATA NPASS/2,  4,   6,   11,   1,   3*0/
C
C     ITYPS DEFINES THE TYPE OF SHELL, 1=S,2=P,3=D
C
      DATA ITYPS/2*1,9*0,   2*1,2*2,7*0,   2*1,2*2,2*3,5*0,
     *           11*1,  1*1,10*0,  11*0,  11*0,   11*0/
C
C     NGS TELLS THE NUMBER OF GAUSSIANS IN EACH SHELL
C
      DATA NGS/1,1,9*0,    1,1,1,1,7*0,  1,1,1,1,1,1,5*0,
     *         1,1,1,1,1,1,1,1,1,1,1,   1,10*0,   11*0,
     *         11*0,   11*0/
C
C     --- PROVIDE "BRAINS" FOR DOUBLE ZETA CONTRACTIONS ---
C
      KIND = KINDS(NN)
      MXPASS = NPASS(KIND)
C
      DO 100 I=1,MXPASS
         ITYPES(I) = ITYPS(I,KIND)
         NGBFS(I) = NGS(I,KIND)
  100 CONTINUE
C
      RETURN
      END
C
C*MODULE NEOBAS  *DECK NDZV01
      SUBROUTINE NDZV01(E,S,P,D,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E(20),S(20),P(20),D(20)
C
      GO TO (100,120,140),N
C
C     ----- HYDROGEN NUCLEUS (2S/2S) -----
C
  100 CONTINUE
      E(1) = 2.9292661D+01
      S(1) = 1.0000000D+00
      E(2) = 2.0175831D+01
      S(2) = 1.0000000D+00
CSPW  E(1) = 2.1100926D+01
CSPW  S(1) = 1.0000000D+00
CSPW  E(2) = 2.0175831D+01
CSPW  S(2) = 1.0000000D+00
      RETURN
C
C     ----- HYDROGEN NUCLEUS (2S6P/2S6P) -----
C
  120 CONTINUE
      E(1) = 2.9292661D+01
      S(1) = 1.0000000D+00
      E(2) = 2.0175831D+01
      S(2) = 1.0000000D+00
      E(3) = 2.6166875D+01
      P(3) = 1.0000000D+00
      E(4) = 1.8217659D+01
      P(4) = 1.0000000D+00
      RETURN
C
C     ----- HYDROGEN NUCLEUS (2S6P12D/2S6P12D) -----
C
  140 CONTINUE
      E(1) = 2.9292661D+01
      S(1) = 1.0000000D+00
      E(2) = 2.0175831D+01
      S(2) = 1.0000000D+00
      E(3) = 2.6166875D+01
      P(3) = 1.0000000D+00
      E(4) = 1.8217659D+01
      P(4) = 1.0000000D+00
      E(5) = 2.3706622D+01
      D(5) = 1.0000000D+00
      E(6) = 1.6730362D+01
      D(6) = 1.0000000D+00
      RETURN
C
      END
C
C*MODULE NEOBAS  *DECK TZNB
      SUBROUTINE TZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,NUCTYP,
     *                IAT,NTSHEL,NNUCSH,NQN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION CSINP(*),CPINP(*),CDINP(*),
     *          INTYP(*),NANGM(*),NBFS(*),MINF(*),MAXF(*),NS(*)
      DIMENSION EEX(31),CCS(31),CCP(31),CCD(31),
     *          ITYPES(11),NGBFS(11)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      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
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, PT5=0.5D+00, PT75=0.75D+00)
      PARAMETER (PI32=5.56832799683170D+00, TOL=1.0D-10)
C
C     ----- SET UP TRIPLE ZETA QUALITY BASIS SET -----
C
      DO 100 I = 1,31
      EEX(I) = ZERO
      CCS(I) = ZERO
      CCP(I) = ZERO
  100 CCD(I) = ZERO
C
C     ----- PROTON -----
C
      IF (NUCTYP.EQ.1) THEN
C
         CALL NTZV01(EEX,CCS,CCP,CCD,NN)
C
C     ----- DEUTERIUM NUCLEUS -----
C
      ELSE IF (NUCTYP.EQ.2) THEN
C
         CALL ABRT
C
C     ----- TRITIUM NUCLEUS -----
C
      ELSE IF (NUCTYP.EQ.3) THEN
C
         CALL ABRT
C
      ELSE
         CALL ABRT
      END IF
C
C     ----- NOW SET UP EACH SHELL -----
C
      CALL NTZSHL(NN,MXPASS,ITYPES,NGBFS)
C
      IGX = 0
      IPASS = 0
  310 CONTINUE
      IPASS = IPASS + 1
      IF(IPASS.GT.MXPASS) RETURN
C                         ******
      IGAUSS = NGBFS(IPASS)
      ITYP = ITYPES(IPASS)
      NNUCSH = NNUCSH+1
      NTSHEL = NTSHEL+1
      IF(NTSHEL.GT.MXSH) THEN
         IERR1=1
         RETURN
      END IF
      NS(NQN) = NS(NQN)+1
      KMIN(NTSHEL) = MINF(ITYP)
      KMAX(NTSHEL) = MAXF(ITYP)
      KSTART(NTSHEL) = NGAUSS+1
      KATOM(NTSHEL) = IAT
      KTYPE(NTSHEL) = NANGM(ITYP)
      INTYP(NTSHEL) = ITYP
      KNG(NTSHEL) = IGAUSS
      KLOC(NTSHEL) = LOC+1
      NGAUSS = NGAUSS+IGAUSS
      IF(NGAUSS.GT.MXGTOT) THEN
         IERR2=1
         RETURN
      END IF
      LOC = LOC+NBFS(ITYP)
      K1 = KSTART(NTSHEL)
      K2 = K1+KNG(NTSHEL)-1
      DO 720 I = 1,IGAUSS
         K = K1+I-1
         EX(K)    = EEX(IGX+I)
         CSINP(K) = CCS(IGX+I)
         CPINP(K) = CCP(IGX+I)
         CDINP(K) = CCD(IGX+I)
         CS(K) = CSINP(K)
         CP(K) = CPINP(K)
         CD(K) = CDINP(K)
  720 CONTINUE
      IGX = IGX + IGAUSS
C
C     ----- ALWAYS UNNORMALIZE PRIMITIVES -----
C
      DO 740 K = K1,K2
         EE = EX(K)+EX(K)
         FACS = PI32/(EE*SQRT(EE))
         FACP = PT5*FACS/EE
         FACD = PT75*FACS/(EE*EE)
         CS(K) = CS(K)/SQRT(FACS)
         CP(K) = CP(K)/SQRT(FACP)
         CD(K) = CD(K)/SQRT(FACD)
  740 CONTINUE
C
C     ----- IF(NORMF.EQ.0) NORMALIZE BASIS FUNCTIONS. -----
C
      IF (NORMF .EQ. 1) GO TO 820
      FACS = ZERO
      FACP = ZERO
      FACD = ZERO
      DO 780 IG = K1,K2
         DO 770 JG = K1,IG
            EE = EX(IG)+EX(JG)
            FAC = EE*SQRT(EE)
            DUMS = CS(IG)*CS(JG)/FAC
            DUMP = PT5*CP(IG)*CP(JG)/(EE*FAC)
            DUMD = PT75*CD(IG)*CD(JG)/(EE*EE*FAC)
            IF (IG .NE. JG) THEN
               DUMS = DUMS+DUMS
               DUMP = DUMP+DUMP
               DUMD = DUMD+DUMD
            END IF
            FACS = FACS+DUMS
            FACP = FACP+DUMP
            FACD = FACD+DUMD
  770    CONTINUE
  780 CONTINUE
      IF (FACS .GT. TOL) FACS = ONE/SQRT(FACS*PI32)
      IF (FACP .GT. TOL) FACP = ONE/SQRT(FACP*PI32)
      IF (FACD .GT. TOL) FACD = ONE/SQRT(FACD*PI32)
      DO 800 IG = K1,K2
         CS(IG) = CS(IG) * FACS
         CP(IG) = CP(IG) * FACP
         CD(IG) = CD(IG) * FACD
  800 CONTINUE
C
C        AND BRANCH BACK FOR THE NEXT SHELL
C
  820 CONTINUE
      GO TO 310
C
      END
C
C*MODULE NEOBAS  *DECK NTZV01
      SUBROUTINE NTZV01(E,S,P,D,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E(20),S(20),P(20),D(20)
C
      GO TO (100,120,140),N
C
C     ----- HYDROGEN NUCLEUS (9S/3S) -----
C
  100 CONTINUE
      E(1) = 6.600000D+01
      S(1) = 1.000000D+00
      E(2) = 1.872000D+01
      S(2) = 1.000000D+00
      E(3) = 5.970000D+00
      S(3) = 1.000000D+00
      RETURN
C
C     ----- HYDROGEN NUCLEUS (3S9P/3S9P) -----
C
  120 CONTINUE
      E(1) = 6.600000D+01
      S(1) = 1.000000D+00
      E(2) = 1.872000D+01
      S(2) = 1.000000D+00
      E(3) = 7.970000D+00
      S(3) = 1.000000D+00
      E(4) = 6.600000D+01
      P(4) = 1.000000D+00
      E(5) = 1.872000D+01
      P(5) = 1.000000D+00
      E(6) = 7.970000D+00
      P(6) = 1.000000D+00
      RETURN
C
C     ----- HYDROGEN NUCLEUS (3S9P18D/3S9P18D) -----
C
  140 CONTINUE
C
      E(1) = 6.600000D+01
      S(1) = 1.000000D+00
      E(2) = 1.872000D+01
      S(2) = 1.000000D+00
      E(3) = 7.970000D+00
      S(3) = 1.000000D+00
      E(4) = 6.600000D+01
      P(4) = 1.000000D+00
      E(5) = 1.872000D+01
      P(5) = 1.000000D+00
      E(6) = 7.970000D+00
      P(6) = 1.000000D+00
      E(7) = 6.600000D+01
      D(7) = 1.000000D+00
      E(8) = 1.872000D+01
      D(8) = 1.000000D+00
      E(9) = 4.970000D+00
      D(9) = 1.000000D+00
      RETURN
      END
C
C*MODULE NEOBAS  *DECK NTZSHL
      SUBROUTINE NTZSHL(NN,MXPASS,ITYPES,NGBFS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION ITYPES(11),NGBFS(11)
      DIMENSION KINDS(36),NPASS(8),ITYPS(11,8),NGS(11,8)
C
C     KINDS ARE S, SP, SPD
C     NPASS TELLS HOW MANY SHELLS ARE IN EACH BASIS
C
      DATA KINDS/1,  2,   3, 33*0/
      DATA NPASS/3,  6,   9,  5*0/
C
C     ITYPS DEFINES THE TYPE OF SHELL, 1=S,2=P,3=D
C
      DATA ITYPS/3*1,8*0,   3*1,3*2,5*0,   3*1,3*2,3*3,2*0,
     *           11*0,  11*0,  11*0,  11*0,   11*0/
C
C     NGS TELLS THE NUMBER OF GAUSSIANS IN EACH SHELL
C
      DATA NGS/1,1,1,8*0,    1,1,1,1,1,1,5*0,  1,1,1,1,1,1,1,1,1,2*0,
     *         11*0,   11*0,   11*0,
     *         11*0,   11*0/
C
C     --- PROVIDE "BRAINS" FOR TRIPLE ZETA CONTRACTIONS ---
C
      KIND = KINDS(NN)
      MXPASS = NPASS(KIND)
C
      DO 100 I=1,MXPASS
         ITYPES(I) = ITYPS(I,KIND)
         NGBFS(I) = NGS(I,KIND)
  100 CONTINUE
C
      RETURN
      END
C
C*MODULE NEOBAS  *DECK NUBLAB
      SUBROUTINE NUBLAB(IARG,ATMSYM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*4 ATMSYM(*)
      CHARACTER*4 BFNAM(35)
      CHARACTER*4 ATMLAB(106)
      CHARACTER*10 BFL
      CHARACTER*8 QNUN,QNN
      CHARACTER*8 LABEL
      CHARACTER*10 PBFLAB
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGTOT=20000, MXAO=8192)
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      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 /NUCMON/ QNUN(20),QNN(20),PBFLAB(MXAO)
C
      DATA BFNAM/'  S ','  X ','  Y ','  Z ',
     *           ' XX ',' YY ',' ZZ ',' XY ',' XZ ',' YZ ',
     *           ' XXX',' YYY',' ZZZ',' XXY',' XXZ',
     *           ' YYX',' YYZ',' ZZX',' ZZY',' XYZ',
     *           'XXXX','YYYY','ZZZZ','XXXY','XXXZ',
     *           'YYYX','YYYZ','ZZZX','ZZZY','XXYY',
     *           'XXZZ','YYZZ','XXYZ','YYXZ','ZZXY'/
C
C                           104 TRUE ELEMENTS, AND 2 DUMMIES
C
      DATA ATMLAB/'H   ','HE  ','LI  ','BE  ','B   ','C   ',
     *            'N   ','O   ','F   ','NE  ','NA  ','MG  ',
     *            'AL  ','SI  ','P   ','S   ','CL  ','AR  ',
     *            'K   ','CA  ','SC  ','TI  ','V   ','CR  ',
     *            'MN  ','FE  ','CO  ','NI  ','CU  ','ZN  ',
     *            'GA  ','GE  ','AS  ','SE  ','BR  ','KR  ',
     *            'RB  ','SR  ','Y   ','ZR  ','NB  ','MO  ',
     *            'TC  ','RU  ','RH  ','PD  ','AG  ','CD  ',
     *            'IN  ','SN  ','SB  ','TE  ','I   ','XE  ',
     *            'CS  ','BA  ','LA  ','CE  ','PR  ','ND  ',
     *            'PM  ','SM  ','EU  ','GD  ','TB  ','DY  ',
     *            'HO  ','ER  ','TM  ','YB  ','LU  ','HF  ',
     *            'TA  ','W   ','RE  ','OS  ','IR  ','PT  ',
     *            'AU  ','HG  ','TL  ','PB  ','BI  ','PO  ',
     *            'AT  ','RN  ','FR  ','RA  ','AC  ','TH  ',
     *            'PA  ','U   ','NP  ','PU  ','AM  ','CM  ',
     *            'BK  ','CF  ','ES  ','FM  ','MD  ','NO  ',
     *            'LR  ','RF  ','X   ','BQ  '/
C
      IF(IARG.EQ.1) GO TO 100
      IF(IARG.EQ.2) GO TO 200
      IF(IARG.EQ.3) GO TO 300
      CALL ABRT
C
C         ----- FILL BASIS FUNCTION SYMBOL TABLE -----
C
  100 CONTINUE
      N = 0
      DO 120 II = 1,NTSHEL
         IAT = KATOM(II)
         LABEL = QNN(1)
         MINI = KMIN(II)
         MAXI = KMAX(II)
         DO 110 I = MINI,MAXI
            N = N+1
            WRITE(UNIT=BFL,FMT='(A4,I2,A4)')
     *            LABEL,MOD(IAT,100),BFNAM(I)
            IF(II.GT.NSHELL) READ (UNIT=BFL,FMT='(A10)') PBFLAB(N)
  110    CONTINUE
  120 CONTINUE
      RETURN
C
C        ----- FILL ATOM SYMBOL TABLE -----
C
  200 CONTINUE
      DO 210 IAT=1,NAT
         NUCZ = INT(ZAN(IAT)) + IZCORE(IAT)
         ATMSYM(IAT) = ATMLAB(105)
         IF(NUCZ.GE.1  .AND.  NUCZ.LE.104) ATMSYM(IAT) = ATMLAB(NUCZ)
  210 CONTINUE
      RETURN
C
C        ----- RETURN THE PERIODIC TABLE -----
C
  300 CONTINUE
      DO 310 IAT=1,106
         ATMSYM(IAT) = ATMLAB(IAT)
  310 CONTINUE
      RETURN
      END
C
C*MODULE NEOBAS  *DECK PREVNU
      SUBROUTINE PREVNU(V,E,M,LDV,ISTART,IEND)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
c
      parameter (MXAO=8192)
C
      CHARACTER*8 QNUN,QNN
      CHARACTER*10 PBFLAB
C
      DIMENSION V(LDV,M),E(M)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NUCMON/ QNUN(20),QNN(20),PBFLAB(MXAO)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- PRINT OUT EIGENDATA (VECTORS AND VALUES) -----
C     THE ROWS WILL BE LABELED WITH THE BASIS FUNCTION TAGS.
C     -V- IS N X M, WITH TRUE LEADING DIMENSION -LDV-
C
      IF (MASWRK) THEN
      MAX = 5
      IF (NPRINT .EQ. 6) MAX = 10
      IMAX = 0
C
  100 IMIN = IMAX+1
      IMAX = IMAX+MAX
      IF (IMAX .GT. M) IMAX = M
      WRITE (IW,9008)
      WRITE (IW,9028) (I,I = IMIN,IMAX)
      WRITE (IW,9008)
      WRITE (IW,9068) (E(I),I = IMIN,IMAX)
      WRITE (IW,9008)
      J = 0
      DO 120 IDX = ISTART,IEND
         J = J + 1
         WRITE (IW,9048) J,PBFLAB(IDX),(V(J,I),I = IMIN,IMAX)
  120 CONTINUE
      IF (IMAX .LT. M) GO TO 100
      ENDIF
      RETURN
C
 9008 FORMAT(1X)
 9028 FORMAT(17X,10(4X,I4,3X))
 9048 FORMAT(I5,2X,A10,10F11.6)
 9068 FORMAT(17X,10F11.4)
      END
C
C*MODULE NEOBAS  *DECK SZNB
      SUBROUTINE SZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,
     *                NUCTYP,IAT,NTSHEL,NNUCSH,NQN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION CSINP(*),CPINP(*),CDINP(*),
     *          INTYP(*),NANGM(*),NBFS(*),MINF(*),MAXF(*),NS(*)
      DIMENSION EEX(30),CCS(30),CCP(30),CCD(30),
     *          ITYPES(11),NGBFS(11)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      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
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, PT5=0.5D+00, PT75=0.75D+00)
      PARAMETER (PI32=5.56832799683170D+00, TOL=1.0D-10)
C
C     ----- SET UP THE "DOUBLE ZETA" NUCLEAR BASIS SET -----
C
      DO 100 I = 1,30
         EEX(I) = ZERO
         CCS(I) = ZERO
         CCP(I) = ZERO
         CCD(I) = ZERO
  100 CONTINUE
C
C     ----- PROTON -----
C
      IF (NUCTYP.EQ.1) THEN
C
         CALL NSZV01(EEX,CCS,CCP,CCD,NN)
C
C     ----- DEUTERIUM NUCLEUS -----
C
      ELSE IF (NUCTYP.EQ.2) THEN
C
         CALL ABRT
C
C     ----- TRITIUM NUCLEUS -----
C
      ELSE IF (NUCTYP.EQ.3) THEN
C
         CALL ABRT
C
      ELSE
         CALL ABRT
      END IF
C
C     ----- NOW SET UP EACH SHELL -----
C
      CALL NSZSHL(NN,MXPASS,ITYPES,NGBFS)
C
      IGX = 0
      IPASS = 0
  310 CONTINUE
      IPASS = IPASS + 1
      IF(IPASS.GT.MXPASS) RETURN
C                         ******
      IGAUSS = NGBFS(IPASS)
      ITYP = ITYPES(IPASS)
C
      NNUCSH = NNUCSH+1
      NTSHEL = NTSHEL+1
      IF(NTSHEL.GT.MXSH) THEN
         IERR1=1
         RETURN
      END IF
      NS(NQN) = NS(NQN)+1
      KMIN(NTSHEL) = MINF(ITYP)
      KMAX(NTSHEL) = MAXF(ITYP)
      KSTART(NTSHEL) = NGAUSS+1
      KATOM(NTSHEL) = IAT
      KTYPE(NTSHEL) = NANGM(ITYP)
      INTYP(NTSHEL) = ITYP
      KNG(NTSHEL) = IGAUSS
      KLOC(NTSHEL) = LOC+1
C
      NGAUSS = NGAUSS+IGAUSS
      IF(NGAUSS.GT.MXGTOT) THEN
         IERR2=1
         RETURN
      END IF
      LOC = LOC+NBFS(ITYP)
      K1 = KSTART(NTSHEL)
      K2 = K1+KNG(NTSHEL)-1
      DO 720 I = 1,IGAUSS
         K = K1+I-1
         EX(K)    = EEX(IGX+I)
         CSINP(K) = CCS(IGX+I)
         CPINP(K) = CCP(IGX+I)
         CDINP(K) = CCD(IGX+I)
         CS(K) = CSINP(K)
         CP(K) = CPINP(K)
         CD(K) = CDINP(K)
  720 CONTINUE
      IGX = IGX + IGAUSS
C
C     ----- ALWAYS UNNORMALIZE PRIMITIVES -----
C
      DO 740 K = K1,K2
         EE = EX(K)+EX(K)
         FACS = PI32/(EE*SQRT(EE))
         FACP = PT5*FACS/EE
         FACD = PT75*FACS/(EE*EE)
         CS(K) = CS(K)/SQRT(FACS)
         CP(K) = CP(K)/SQRT(FACP)
         CD(K) = CD(K)/SQRT(FACD)
  740 CONTINUE
C
C     ----- IF(NORMF.EQ.0) NORMALIZE BASIS FUNCTIONS. -----
C
      IF (NORMF .EQ. 1) GO TO 820
      FACS = ZERO
      FACP = ZERO
      FACD = ZERO
      DO 780 IG = K1,K2
         DO 770 JG = K1,IG
            EE = EX(IG)+EX(JG)
            FAC = EE*SQRT(EE)
            DUMS = CS(IG)*CS(JG)/FAC
            DUMP = PT5*CP(IG)*CP(JG)/(EE*FAC)
            DUMD = PT75*CD(IG)*CD(JG)/(EE*EE*FAC)
            IF (IG .NE. JG) THEN
               DUMS = DUMS+DUMS
               DUMP = DUMP+DUMP
               DUMD = DUMD+DUMD
            END IF
            FACS = FACS+DUMS
            FACP = FACP+DUMP
            FACD = FACD+DUMD
  770    CONTINUE
  780 CONTINUE
      IF (FACS .GT. TOL) FACS = ONE/SQRT(FACS*PI32)
      IF (FACP .GT. TOL) FACP = ONE/SQRT(FACP*PI32)
      IF (FACD .GT. TOL) FACD = ONE/SQRT(FACD*PI32)
      DO 800 IG = K1,K2
         CS(IG) = CS(IG) * FACS
         CP(IG) = CP(IG) * FACP
         CD(IG) = CD(IG) * FACD
  800 CONTINUE
C
C        AND BRANCH BACK FOR THE NEXT SHELL
C
  820 CONTINUE
      GO TO 310
      END
C
C*MODULE NEOBAS  *DECK NSZSHL
      SUBROUTINE NSZSHL(NN,MXPASS,ITYPES,NGBFS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION ITYPES(11),NGBFS(11)
      DIMENSION KINDS(36),NPASS(8),ITYPS(11,8),NGS(11,8)
C
C     KINDS ARE S, SP, SPD
C     NPASS TELLS HOW MANY SHELLS ARE IN EACH BASIS
C
      DATA KINDS/1,  2,   3,   4,  5,  31*0/
      DATA NPASS/1,  2,   3,   11,   1,   3*0/
C
C     ITYPS DEFINES THE TYPE OF SHELL, 1=S,2=P,3=D
C
      DATA ITYPS/1*1,10*0,   1*1,1*2,9*0,   1*1,1*2,1*3,8*0,
     *           11*1,  1*1,10*0,  11*0,  11*0,   11*0/
C
C     NGS TELLS THE NUMBER OF GAUSSIANS IN EACH SHELL
C
      DATA NGS/1,1,9*0,    1,1,1,1,7*0,  1,1,1,1,1,1,5*0,
     *         1,1,1,1,1,1,1,1,1,1,1,   1,10*0,   11*0,
     *         11*0,   11*0/
C
C     --- PROVIDE "BRAINS" FOR DOUBLE ZETA CONTRACTIONS ---
C
      KIND = KINDS(NN)
      MXPASS = NPASS(KIND)
C
      DO 100 I=1,MXPASS
         ITYPES(I) = ITYPS(I,KIND)
         NGBFS(I) = NGS(I,KIND)
  100 CONTINUE
C
      RETURN
      END
C
C*MODULE NEOBAS  *DECK NSZV01
      SUBROUTINE NSZV01(E,S,P,D,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E(20),S(20),P(20),D(20)
C
      GO TO (100,120,140),N
C
C     ----- HYDROGEN NUCLEUS (2S/2S) -----
C
  100 CONTINUE
      E(1) = 2.9292661D+01
      S(1) = 1.0000000D+00
      RETURN
C
C     ----- HYDROGEN NUCLEUS (2S6P/2S6P) -----
C
  120 CONTINUE
      E(1) = 2.9292661D+01
      S(1) = 1.0000000D+00
      E(2) = 2.6166875D+01
      P(2) = 1.0000000D+00
      RETURN
C
C     ----- HYDROGEN NUCLEUS (2S6P12D/2S6P12D) -----
C
  140 CONTINUE
      E(1) = 2.9292661D+01
      S(1) = 1.0000000D+00
      E(2) = 2.6166875D+01
      P(2) = 1.0000000D+00
      E(3) = 2.3706622D+01
      D(3) = 1.0000000D+00
      RETURN
C
      END
C
C*MODULE NEOBAS  *DECK QZNB
      SUBROUTINE QZNB(NN,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,
     *                NUCTYP,IAT,NTSHEL,NNUCSH,NQN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION CSINP(*),CPINP(*),CDINP(*),
     *          INTYP(*),NANGM(*),NBFS(*),MINF(*),MAXF(*),NS(*)
      DIMENSION EEX(30),CCS(30),CCP(30),CCD(30),
     *          ITYPES(15),NGBFS(15)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      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
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, PT5=0.5D+00, PT75=0.75D+00)
      PARAMETER (PI32=5.56832799683170D+00, TOL=1.0D-10)
C
C     ----- SET UP THE "QUADRUPLE ZETA" NUCLEAR BASIS SET -----
C
      DO 100 I = 1,30
         EEX(I) = ZERO
         CCS(I) = ZERO
         CCP(I) = ZERO
         CCD(I) = ZERO
  100 CONTINUE
C
C     ----- PROTON -----
C
      IF (NUCTYP.EQ.1) THEN
C
         CALL NQZV01(EEX,CCS,CCP,CCD,NN)
C
C     ----- DEUTERIUM NUCLEUS -----
C
      ELSE IF (NUCTYP.EQ.2) THEN
C
         CALL ABRT
C
C     ----- TRITIUM NUCLEUS -----
C
      ELSE IF (NUCTYP.EQ.3) THEN
C
         CALL ABRT
C
      ELSE
         CALL ABRT
      END IF
C
C     ----- NOW SET UP EACH SHELL -----
C
      CALL NQZSHL(NN,MXPASS,ITYPES,NGBFS)
C
      IGX = 0
      IPASS = 0
  310 CONTINUE
      IPASS = IPASS + 1
      IF(IPASS.GT.MXPASS) RETURN
C                         ******
      IGAUSS = NGBFS(IPASS)
      ITYP = ITYPES(IPASS)
C
      NNUCSH = NNUCSH+1
      NTSHEL = NTSHEL+1
      IF(NTSHEL.GT.MXSH) THEN
         IERR1=1
         RETURN
      END IF
      NS(NQN) = NS(NQN)+1
      KMIN(NTSHEL) = MINF(ITYP)
      KMAX(NTSHEL) = MAXF(ITYP)
      KSTART(NTSHEL) = NGAUSS+1
      KATOM(NTSHEL) = IAT
      KTYPE(NTSHEL) = NANGM(ITYP)
      INTYP(NTSHEL) = ITYP
      KNG(NTSHEL) = IGAUSS
      KLOC(NTSHEL) = LOC+1
C
      NGAUSS = NGAUSS+IGAUSS
      IF(NGAUSS.GT.MXGTOT) THEN
         IERR2=1
         RETURN
      END IF
      LOC = LOC+NBFS(ITYP)
      K1 = KSTART(NTSHEL)
      K2 = K1+KNG(NTSHEL)-1
      DO 720 I = 1,IGAUSS
         K = K1+I-1
         EX(K)    = EEX(IGX+I)
         CSINP(K) = CCS(IGX+I)
         CPINP(K) = CCP(IGX+I)
         CDINP(K) = CCD(IGX+I)
         CS(K) = CSINP(K)
         CP(K) = CPINP(K)
         CD(K) = CDINP(K)
  720 CONTINUE
      IGX = IGX + IGAUSS
C
C     ----- ALWAYS UNNORMALIZE PRIMITIVES -----
C
      DO 740 K = K1,K2
         EE = EX(K)+EX(K)
         FACS = PI32/(EE*SQRT(EE))
         FACP = PT5*FACS/EE
         FACD = PT75*FACS/(EE*EE)
         CS(K) = CS(K)/SQRT(FACS)
         CP(K) = CP(K)/SQRT(FACP)
         CD(K) = CD(K)/SQRT(FACD)
  740 CONTINUE
C
C     ----- IF(NORMF.EQ.0) NORMALIZE BASIS FUNCTIONS. -----
C
      IF (NORMF .EQ. 1) GO TO 820
      FACS = ZERO
      FACP = ZERO
      FACD = ZERO
      DO 780 IG = K1,K2
         DO 770 JG = K1,IG
            EE = EX(IG)+EX(JG)
            FAC = EE*SQRT(EE)
            DUMS = CS(IG)*CS(JG)/FAC
            DUMP = PT5*CP(IG)*CP(JG)/(EE*FAC)
            DUMD = PT75*CD(IG)*CD(JG)/(EE*EE*FAC)
            IF (IG .NE. JG) THEN
               DUMS = DUMS+DUMS
               DUMP = DUMP+DUMP
               DUMD = DUMD+DUMD
            END IF
            FACS = FACS+DUMS
            FACP = FACP+DUMP
            FACD = FACD+DUMD
  770    CONTINUE
  780 CONTINUE
      IF (FACS .GT. TOL) FACS = ONE/SQRT(FACS*PI32)
      IF (FACP .GT. TOL) FACP = ONE/SQRT(FACP*PI32)
      IF (FACD .GT. TOL) FACD = ONE/SQRT(FACD*PI32)
      DO 800 IG = K1,K2
         CS(IG) = CS(IG) * FACS
         CP(IG) = CP(IG) * FACP
         CD(IG) = CD(IG) * FACD
  800 CONTINUE
C
C        AND BRANCH BACK FOR THE NEXT SHELL
C
  820 CONTINUE
      GO TO 310
      END
C
C*MODULE NEOBAS  *DECK NQZSHL
      SUBROUTINE NQZSHL(NN,MXPASS,ITYPES,NGBFS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION ITYPES(15),NGBFS(15)
      DIMENSION KINDS(36),NPASS(8),ITYPS(15,8),NGS(15,8)
C
C     KINDS ARE S, SP, SPD
C     NPASS TELLS HOW MANY SHELLS ARE IN EACH BASIS
C
      DATA KINDS/1,  2,   3,   33*0/
      DATA NPASS/4,  8,  10,    5*0/
C
C     ITYPS DEFINES THE TYPE OF SHELL, 1=S,2=P,3=D
C
      DATA ITYPS/4*1,11*0,   4*1,4*2,7*0,   4*1,4*2,2*3,5*0,
     *           15*0,  15*0,  15*0,  15*0,   15*0/
C
C     NGS TELLS THE NUMBER OF GAUSSIANS IN EACH SHELL
C
      DATA NGS/1,1,1,1,11*0,    1,1,1,1,1,1,1,1,7*0,
     *         1,1,1,1,1,1,1,1,1,1,5*0,    15*0,   15*0,   15*0,
     *         15*0,   15*0/
C
C     --- PROVIDE "BRAINS" FOR QUADRUPLE ZETA CONTRACTIONS ---
C
      KIND = KINDS(NN)
      MXPASS = NPASS(KIND)
C
      DO 100 I=1,MXPASS
         ITYPES(I) = ITYPS(I,KIND)
         NGBFS(I) = NGS(I,KIND)
  100 CONTINUE
C
      RETURN
      END
C
C*MODULE NEOBAS  *DECK NQZV01
      SUBROUTINE NQZV01(E,S,P,D,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E(20),S(20),P(20),D(20)
C
      GO TO (100,120,140),N
C
C     ----- HYDROGEN NUCLEUS (4S/4S) -----
C
  100 CONTINUE
      E(1) = 3.463529D+01
      S(1) = 1.000000D+00
      E(2) = 2.353722D+01
      S(2) = 1.000000D+00
      E(3) = 1.615827D+01
      S(3) = 1.000000D+00
      E(4) = 2.880598D+00
      S(4) = 1.000000D+00
      RETURN
C
C     ----- HYDROGEN NUCLEUS (4S12P/4S12P) -----
C
  120 CONTINUE
      E(1) = 3.463529D+01
      S(1) = 1.000000D+00
      E(2) = 2.353722D+01
      S(2) = 1.000000D+00
      E(3) = 1.615827D+01
      S(3) = 1.000000D+00
      E(4) = 2.880598D+00
      S(4) = 1.000000D+00
      E(5) = 3.360130D+01
      P(5) = 1.000000D+00
      E(6) = 2.291236D+01
      P(6) = 1.000000D+00
      E(7) = 1.580526D+01
      P(7) = 1.000000D+00
      E(8) = 6.897948D+00
      P(8) = 1.000000D+00
      RETURN
C
C     ----- HYDROGEN NUCLEUS (4S12P12D/4S12P12D) -----
C
  140 CONTINUE
      E(1) = 3.463529D+01
      S(1) = 1.000000D+00
      E(2) = 2.353722D+01
      S(2) = 1.000000D+00
      E(3) = 1.615827D+01
      S(3) = 1.000000D+00
      E(4) = 2.880598D+00
      S(4) = 1.000000D+00
      E(5) = 3.360130D+01
      P(5) = 1.000000D+00
      E(6) = 2.291236D+01
      P(6) = 1.000000D+00
      E(7) = 1.580526D+01
      P(7) = 1.000000D+00
      E(8) = 6.897948D+00
      P(8) = 1.000000D+00
      E(9) = 2.444627D+01
      D(9) = 1.000000D+00
      E(10) = 1.712008D+01
      D(10) = 1.000000D+00
      RETURN
C
      END
C
C*MODULE NEOBAS  *DECK READNB
      SUBROUTINE READNB(IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *                INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,
     *                IAT,NTSHEL,NNUCSH,NQN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
CSPW  CHARACTER*(*) NUCBAS
CSPW  CHARACTER*256 NUBPTH
C
      DIMENSION CSINP(*),CPINP(*),CDINP(*),
     *          INTYP(*),NANGM(*),NBFS(*),MINF(*),MAXF(*),NS(*)
      DIMENSION EEX(30),CCS(30),CCP(30),CCD(30),
     *          ITYPES(12),NGBFS(12)
C
      PARAMETER (MXSH=5000, MXGTOT=20000)
C
      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
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, PT5=0.5D+00, PT75=0.75D+00)
      PARAMETER (PI32=5.56832799683170D+00, TOL=1.0D-10)
C
      DO 100 I = 1,30
         EEX(I) = ZERO
         CCS(I) = ZERO
         CCP(I) = ZERO
         CCD(I) = ZERO
  100 CONTINUE
C
C     ----- PROTON -----
C
      IBF=3
      CALL SEQOPN(IBF,'NUCBAS', 'OLD',.TRUE., 'FORMATTED')
      READ(IBF,*) MXPASS
      DO 200 I=1,MXPASS
         NGBFS(I) = 1
         READ(IBF,*) ITYPES(I), EEX(I)
         IF (ITYPES(I).EQ.1) CCS(I)= ONE
         IF (ITYPES(I).EQ.2) CCP(I)= ONE
         IF (ITYPES(I).EQ.3) CCD(I)= ONE
  200 CONTINUE
C
      IGX = 0
      IPASS = 0
  310 CONTINUE
      IPASS = IPASS + 1
      IF(IPASS.GT.MXPASS) RETURN
C                         ******
      IGAUSS = NGBFS(IPASS)
      ITYP = ITYPES(IPASS)
C
      NNUCSH = NNUCSH+1
      NTSHEL = NTSHEL+1
      IF(NTSHEL.GT.MXSH) THEN
         IERR1=1
         RETURN
      END IF
      NS(NQN) = NS(NQN)+1
      KMIN(NTSHEL) = MINF(ITYP)
      KMAX(NTSHEL) = MAXF(ITYP)
      KSTART(NTSHEL) = NGAUSS+1
      KATOM(NTSHEL) = IAT
      KTYPE(NTSHEL) = NANGM(ITYP)
      INTYP(NTSHEL) = ITYP
      KNG(NTSHEL) = IGAUSS
      KLOC(NTSHEL) = LOC+1
      NGAUSS = NGAUSS+IGAUSS
      IF(NGAUSS.GT.MXGTOT) THEN
           IERR2=1
           RETURN
      END IF
      LOC = LOC+NBFS(ITYP)
      K1 = KSTART(NTSHEL)
      K2 = K1+KNG(NTSHEL)-1
      DO 720 I = 1,IGAUSS
            K = K1+I-1
             EX(K)    = EEX(IGX+I)
              CSINP(K) = CCS(IGX+I)
              CPINP(K) = CCP(IGX+I)
         CDINP(K) = CCD(IGX+I)
          CS(K) = CSINP(K)
          CP(K) = CPINP(K)
          CD(K) = CDINP(K)
  720 CONTINUE
      IGX = IGX + IGAUSS
C
C     ----- ALWAYS UNNORMALIZE PRIMITIVES -----
C
      DO 740 K = K1,K2
         EE = EX(K)+EX(K)
         FACS = PI32/(EE*SQRT(EE))
         FACP = PT5*FACS/EE
         FACD = PT75*FACS/(EE*EE)
         CS(K) = CS(K)/SQRT(FACS)
         CP(K) = CP(K)/SQRT(FACP)
         CD(K) = CD(K)/SQRT(FACD)
  740 CONTINUE
C
C     ----- IF(NORMF.EQ.0) NORMALIZE BASIS FUNCTIONS. -----
C
      IF (NORMF .EQ. 1) GO TO 820
      FACS = ZERO
      FACP = ZERO
      FACD = ZERO
      DO 780 IG = K1,K2
         DO 770 JG = K1,IG
            EE = EX(IG)+EX(JG)
            FAC = EE*SQRT(EE)
            DUMS = CS(IG)*CS(JG)/FAC
            DUMP = PT5*CP(IG)*CP(JG)/(EE*FAC)
            DUMD = PT75*CD(IG)*CD(JG)/(EE*EE*FAC)
           IF (IG .NE. JG) THEN
              DUMS = DUMS+DUMS
              DUMP = DUMP+DUMP
              DUMD = DUMD+DUMD
           END IF
           FACS = FACS+DUMS
           FACP = FACP+DUMP
           FACD = FACD+DUMD
  770    CONTINUE
  780 CONTINUE
      IF (FACS .GT. TOL) FACS = ONE/SQRT(FACS*PI32)
      IF (FACP .GT. TOL) FACP = ONE/SQRT(FACP*PI32)
      IF (FACD .GT. TOL) FACD = ONE/SQRT(FACD*PI32)
      DO 800 IG = K1,K2
         CS(IG) = CS(IG) * FACS
         CP(IG) = CP(IG) * FACP
         CD(IG) = CD(IG) * FACD
  800 CONTINUE
C
C        AND BRANCH BACK FOR THE NEXT SHELL
C
  820 CONTINUE
      GO TO 310
      END
C
