C  7 DEC 07 - MWS - DEVELOP AUTOMATED BONDING ANALYSIS PATHWAY
C 10 JUL 06 - MWS - FIX MBS COUNTING ESTIMATES
C  8 MAY 06 - MWS - NEW MODULE FOR QUAMBO'S/VIRTUAL VALENCE ORBITALS
C
C*MODULE VVOS    *DECK VVOS
      SUBROUTINE VVOS(VEC,FAO,EIG,WRK1,WRK2,WRK3,SCR,IWRK,
     *                L0,L1,L2,L3,NOCC)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION VEC(L1,L0),FAO(L2),EIG(L1),
     *          WRK1(L3),WRK2(L3),WRK3(L3),SCR(L1,8),IWRK(L1)
      DIMENSION NBFS(4),MINF(4),MAXF(4),NANGM(4)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,OK,DZ
C
      PARAMETER (MXGTOT=20000, MXSH=5000, MXATM=2000)
C
      COMMON /FMCOM / XX(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NEX,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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TT(432),INVT(48),NT
      COMMON /VVOPAR/ IVVOS,NVVOS,BNDDEN
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C         THE PERIODIC TABLE ONLY GOES UP TO SPDF:
      DATA NBFS/ 1, 3, 6,10/
      DATA MINF/ 1, 2, 5,11/
      DATA MAXF/ 1, 4,10,20/
      DATA NANGM/1, 2, 3, 4/
      DATA CHECK/8HCHECK   /
C
C         ----- GENERATE VIRTUAL VALENCE ORBITAL SPACE -----
C       SEE PARTICULARLY SECTION II.F AND V OF THE FIRST PAPER:
C
C     MOLECULE INTRINSIC MINIMAL BASIS SETS. I. EXACT RESOLUTION OF
C     AB INITIO OPTIMIZED MOLECULAR ORBITALS IN TERMS OF DEFORMED
C     MINIMAL-BASIS ORBITALS
C     W.C.LU, C.Z.WANG, M.W.SCHMIDT, L.BYTAUTAS, K.M.HO, K.RUEDENBERG
C     J.CHEM.PHYS. 120, 2629-2637(2004)
C
C     MOLECULE INTRINSIC MINIMAL BASIS SETS. II. BONDING ANALYSIS
C     FOR SI4H6 AND SI2 TO SI10.
C     W.C.LU, C.Z.WANG, M.W.SCHMIDT, L.BYTAUTAS, K.M.HO, K.RUEDENBERG
C     J.CHEM.PHYS. 120, 2638-2651(2004)
C
C     REPRESENTATION OF ELECTRONIC STRUCTURES IN CRYSTALS IN TERMS OF
C     HIGHLY LOCALIZED QUASIATOMIC MINIMAL BASIS ORBITALS.
C     W.C.LU, C.Z.WANG, T.L.CHAN, K.RUEDENBERG, K.M.HO
C     PHYS.REV.B 70, 041101/1-4(2004)
C
C     ON ENTRY, -FAO- CONTAINS THE RHF FOCK OPERATOR IN THE AO BASIS,
C     WHILE OTHER ARRAYS ARE TREATED AS WORK SPACES.
C
C     BE SURE THE RUN IS AN ALL E- BASIS SET AND ONLY H-AR PRESENT,
C     BECAUSE THESE ARE THE ONLY KIND OF ATOMIC MINIMAL BASIS SETS
C     THAT ARE INTERNALLY STORED FOR USE IN THE PROJECTION OF AOS.
C
      OK = .TRUE.
      DO I=1,NAT
         NUCZ = INT(ZAN(I) + 0.001D+00)
         IF(IAN(I).GT.18)   OK=.FALSE.
         IF(NUCZ.NE.IAN(I)) OK=.FALSE.
      ENDDO
      IF(.NOT.OK) THEN
         IF(MASWRK) WRITE(IW,9000)
         RETURN
      END IF
C
      NVALAT = 0
      DO I=1,NAT
         IF(                     IAN(I).LE. 2) NVALAT=NVALAT+1
         IF(IAN(I).GE. 3  .AND.  IAN(I).LE. 4) NVALAT=NVALAT+1
         IF(IAN(I).GE. 5  .AND.  IAN(I).LE.10) NVALAT=NVALAT+4
         IF(IAN(I).GE.11  .AND.  IAN(I).LE.12) NVALAT=NVALAT+1
         IF(IAN(I).GE.13  .AND.  IAN(I).LE.18) NVALAT=NVALAT+4
      ENDDO
      NCORE   = NUMCOR()
      NOCCVAL = NOCC - NCORE
      NVIRVAL = NVALAT - NOCCVAL
      NVIR    = L0 - NOCC
      IF(MASWRK) WRITE(IW,9010) NCORE,NVALAT,NOCCVAL,NVIRVAL,NVIR
C
      NVVOS = NVIRVAL
C
C         ESTIMATE GAUSSIANS AND SHELLS FOR BIGGEST ATOM IN EACH ROW.
      MXGMBS = 0
      MXSHMBS= 0
      DO IAT=1,NAT
         IF(                       IAN(IAT).LE. 2) THEN
            MXGMBS  = MXGMBS  + 8
            MXSHMBS = MXSHMBS + 1
         END IF
         IF(IAN(IAT).GE. 3  .AND.  IAN(IAT).LE.10) THEN
            MXGMBS  = MXGMBS  + 14 + 14 + 7
            MXSHMBS = MXSHMBS + 3
         END IF
         IF(IAN(IAT).GE.11  .AND.  IAN(IAT).LE.18) THEN
            MXGMBS  = MXGMBS  + 18 + 18 + 12 + 18 + 12
            MXSHMBS = MXSHMBS + 5
         END IF
      ENDDO
C
      CALL VALFM(LOADFM)
      LEXMBS   = LOADFM   + 1
      LCSMBS   = LEXMBS   + MXGMBS
      LCPMBS   = LCSMBS   + MXGMBS
      LCSINP   = LCPMBS   + MXGMBS
      LCPINP   = LCSINP   + MXGMBS
      LINTYP   = LCPINP   + MXGMBS
      LKSTRMBS = LINTYP   + MXSHMBS
      LKATMMBS = LKSTRMBS + MXSHMBS
      LKTYPMBS = LKATMMBS + MXSHMBS
      LKNGMBS  = LKTYPMBS + MXSHMBS
      LKLOCMBS = LKNGMBS  + MXSHMBS
      LKMINMBS = LKLOCMBS + MXSHMBS
      LKMAXMBS = LKMINMBS + MXSHMBS
      LNSHLS   = LKMAXMBS + MXSHMBS
      LAST     = LNSHLS   + NAT
      NEED1 = LAST -LOADFM -1
      CALL GETFM(NEED1)
C
C     ---- SET UP NEAR-HF QUALITY MINIMAL ATOMIC BASIS, ON THE FLY ----
C     BESIDES SETTING UP DETAILS IN EXMBS, CSMBS, CPMBS, KXXXMBS, ...,
C     THE DIMENSIONALITY IS SPECIFIED BY L1MBS,NGAUMBS,NSHLMBS
C     NOTE THAT CSINP,CPINP,INTYP,NSHLS ARE STRICTLY WORKING STORAGE.
C     CAUTION, NAT NEEDS TO COUNT 1,2,3,... AS THE BASIS IS CONSTRUCTED.
C
      IERR1   = 0
      IERR2   = 0
      L1MBS   = 0
      NSHLMBS = 0
      NGAUMBS = 0
      DZ = .FALSE.
      NATSV = NAT
C
      DO IAT=1,NAT
         NAT  = IAT
         NUCZ = IAN(IAT)
         CALL ETGTO(NUCZ,DZ,XX(LCSINP),XX(LCPINP),IERR1,IERR2,
     *              L1MBS,XX(LINTYP),NANGM,NBFS,MINF,MAXF,XX(LNSHLS),
     *              XX(LEXMBS),XX(LCSMBS),XX(LCPMBS),XX(LKSTRMBS),
     *              XX(LKATMMBS),XX(LKTYPMBS),XX(LKNGMBS),
     *              XX(LKLOCMBS),XX(LKMINMBS),XX(LKMAXMBS),
     *              NGAUMBS,NSHLMBS,MXGMBS,MXSHMBS)
         IF(MASWRK  .AND.  IERR1.NE.0) WRITE(IW,9020) MXSHMBS
         IF(MASWRK  .AND.  IERR2.NE.0) WRITE(IW,9030) MXGMBS
         IF(IERR1.NE.0  .OR.  IERR2.NE.0) CALL ABRT
      ENDDO
      NAT = NATSV
C
      CALL VALFM(LOADFM)
      LOVLP = LOADFM + 1
      LAST  = LOVLP  + L1*L1MBS
      NEED2 = LAST -LOADFM -1
      CALL GETFM(NEED2)
      IF(EXETYP.EQ.CHECK) GO TO 110
C
C        STEP (I), GENERATION OF THE MATRIX "LITTLE-A-STAR", IN -WRK1-
C
C        THE IMPLEMENTATION HERE USES A FIXED, PREVIOUSLY STORED ATOMIC
C        BASIS SET AS THE CAPITAL-A-STAR ORBITALS, RATHER THAN TAKING
C        THESE FROM AN ATOMIC COMPUTATION USING THE CURRENT MOLECULE'S
C        BASIS SET.
C
C        GET OVERLAP BETWEEN THE AO MBS AND THE BASIS SET OF THE RUN,
C        COOVLP COMPUTES <BASIS2|BASIS1>, WITH SHAPE = OVLP(L1,L1MBS),
C        THEN TRANSFORM FROM AO TO THE SCF MO BASIS (OCC+VIRT).
C
      MODE=1
      DUMMY=ZERO
      IDUMMY=1
      CALL COOVLP(MODE,DUMMY,XX(LOVLP),L1MBS,L1,IDUMMY,NAT,
     *            NGAUMBS,NSHLMBS,XX(LEXMBS),XX(LCSMBS),XX(LCPMBS),
     *            DUMMY,DUMMY,DUMMY,DUMMY,DUMMY,
     *            XX(LKSTRMBS),XX(LKATMMBS),XX(LKTYPMBS),XX(LKNGMBS),
     *            XX(LKLOCMBS),XX(LKMINMBS),XX(LKMAXMBS),
     *            MXGTOT,NSHELL,EX,CS,CP,CD,CF,CG,CH,CI,
     *            KSTART,KATOM,KTYPE,KNG,KLOC,KMIN,KMAX,
     *            C,C)
      CALL DGEMM('T','N',L0,L1MBS,L1,
     *           ONE,VEC,L1,XX(LOVLP),L1,ZERO,WRK1,L0)
C
  110 CONTINUE
      CALL RETFM(NEED2)
      CALL RETFM(NEED1)
      IF(EXETYP.EQ.CHECK) GO TO 310
C
C        STEP (II), FORM CAPITAL-B (IN -WRK2-) AND DIAGONALIZE IT
C                   TO OBTAIN ITS -P- (NVIRVAL) LARGEST VALUES
C
      IJ=0
      DO I=1,NVIR
         DO J=1,I
            IJ=IJ+1
            WRK2(IJ) = DDOT(L1MBS,WRK1(NOCC+I),L0,WRK1(NOCC+J),L0)
         ENDDO
      ENDDO
C         DIAGONALIZATION SHOULD GET EIGENVALUES IN DESCENDING ORDER
      NVIR2 = (NVIR*NVIR+NVIR)/2
      DO I=1,NVIR2
         WRK2(I) = -WRK2(I)
      ENDDO
      CALL GLDIAG(NVIR,NVIR,NVIR,WRK2,SCR,EIG,WRK3,IERR,IWRK)
      DO I=1,NVIR
         EIG(I) = -EIG(I)
      ENDDO
C
C        SHOW THE USER THE BREAK IN EIGENVALUES FOR -B-.
C        PERHAPS SOMEDAY A NUMERICAL TEST CAN BE DEVISED HERE, ONCE
C        SOME EXPERIENCE IN THE SIZES OF EIGENVALUES IS GAINED.
C
      MINPRT=MAX(   1,NVIRVAL-5)
      MAXPRT=MIN(NVIR,NVIRVAL+5)
      IF(MASWRK) THEN
         WRITE(IW,9040)
         WRITE(IW,9050)
         WRITE(IW,9070) (EIG(III),III=MINPRT,NVIRVAL)
         WRITE(IW,9060)
         WRITE(IW,9070) (EIG(III),III=NVIRVAL+1,MAXPRT)
      END IF
C
C        STEP (III), FORM CAPITAL-R, CAPITAL-D, BOTH LITTLE-A-PRIME'S
C
  310 CONTINUE
      CALL VALFM(LOADFM)
      LR    = LOADFM + 1
      LD    = LR     + NVIR*NVIR
      LWRK4 = LD     + L1MBS
      LWRK5 = LWRK4  + NVIR*L1MBS
      LAST  = LWRK5  + L1MBS*L1MBS
      NEED = LAST-LOADFM-1
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 500
C
C          NOTE THAT THE PROJECTOR R FORMED HERE SATISFIES R*R=R
C
      CALL DGEMM('N','T',NVIR,NVIR,NVIRVAL,
     *           ONE,WRK3,NVIR,WRK3,NVIR,ZERO,XX(LR),NVIR)
C
C          COMPUTE 2ND TERM OF D BEFORE THE 1ST TERM,
C          NOTE REUSE OF THE INTERMEDIATE -WRK4- IN THE SECOND LOOP.
C
      CALL DGEMM('N','N',NVIR,L1MBS,NVIR,
     *           ONE,XX(LR),NVIR,WRK1(NOCC+1),L0,ZERO,XX(LWRK4),NVIR)
      CALL DGEMM('T','N',L1MBS,1,NVIR,
     *           ONE,WRK1(NOCC+1),L0,XX(LWRK4),NVIR,ZERO,XX(LD),L1MBS)
      DO J=1,L1MBS
         LW1 = L0*(J-1) + 1
         XX(LD+J-1) = XX(LD+J-1) + DDOT(NOCC,WRK1(LW1),1,WRK1(LW1),1)
      ENDDO
C
C          BUILD LITTLE-A-PRIME ON TOP OF LITTLE-A-STAR, IN -WRK1-
C
      DO J=1,L1MBS
         LW1 =   L0*(J-1) + 1
         LW4 = NVIR*(J-1) + LWRK4
         DFACT = ONE/SQRT(XX(LD+J-1))
         CALL DSCAL(NOCC,DFACT,WRK1(LW1),1)
         CALL DSCAL(NVIR,DFACT,  XX(LW4),1)
         CALL DCOPY(NVIR,XX(LW4),1,WRK1(LW1+NOCC),1)
      ENDDO
C
C        STEP (V), SYMMETRIC ORTHOGONALIZATION ACROSS ALL ATOMS,
C                  THUS PRODUCING THE QUAMBO'S, A''.
C
C          FORM OVERLAP MATRIX OVER THE A-PRIME ORBITALS
C
      CALL DGEMM('T','N',L1MBS,L1MBS,NOCC,
     *           ONE,WRK1,L0,WRK1,L0,ZERO,XX(LWRK5),L1MBS)
      CALL DGEMM('T','N',L1MBS,L1MBS,NVIR,
     *           ONE,WRK1(NOCC+1),L0,WRK1(NOCC+1),L0,
     *           ONE,XX(LWRK5),L1MBS)
C
C          FORM S**(-1/2) MATRIX, IN -WRK3-
C
      CALL CPYSQT(XX(LWRK5),WRK2,L1MBS,1)
      CALL GLDIAG(L1MBS,L1MBS,L1MBS,WRK2,SCR,EIG,WRK3,IERR,IWRK)
      L2MBS = (L1MBS*L1MBS+L1MBS)
      CALL VCLR(WRK2,1,L2MBS)
      II = 0
      DO I=1,L1MBS
         II=II+I
         WRK2(II) = ONE/SQRT(EIG(I))
      ENDDO
      CALL TRPOSQ(WRK3,L1MBS)
      CALL TFTRI(XX(LWRK5),WRK2,WRK3,SCR,L1MBS,L1MBS,L1MBS)
      CALL CPYTSQ(XX(LWRK5),WRK3,L1MBS,1)
C
C          COMPUTE A'' COEFFICIENTS, A'' = LITTLE-A-PRIME * S**(-1/2)
C          AND THEN COPY THEM TO -WRK1- (OVERWRITING LITTLE-A-PRIME).
C
      CALL DGEMM('N','N',L0,L1MBS,L1MBS,
     *           ONE,WRK1,L0,WRK3,L1MBS,ZERO,WRK2,L0)
      CALL DCOPY(L0*L1MBS,WRK2,1,WRK1,1)
C
C          TRANSFORM A'' FROM CANONICAL MO BASIS TO AO SPACE
C          AND THEN COPY THEM TO -WRK1-
C
      CALL DAREAD(IDAF,IODA,VEC,L3,15,0)
      CALL DGEMM('N','N',L1,L1MBS,L0,
     *           ONE,VEC,L1,WRK1,L0,ZERO,WRK2,L1)
      CALL DCOPY(L1*L1MBS,WRK2,1,WRK1,1)
C
C          DUMP THE QUAMBO'S TO THE PUNCH FILE
C
      IF(MASWRK) THEN
         WRITE(IW,9080) L1MBS,NOCC,NVIRVAL
         WRITE(IP,8000)
         WRITE(IP,8010) ' $VEC'
         CALL PUSQL(WRK1,L1MBS,L1,L1)
         WRITE(IP,8010) ' $END'
      END IF
C
C        STEP (VI), FORM CANONICAL VALENCE ORBITALS
C          THE QUAMBO'S ARE IN WRK1, AND WE ARE MISSING THE
C          EXTERNAL SET.  FORTUNATELY, THE -ORTHO- ROUTINE
C          CAN GENERATE THESE BY ORTHOGONALIZING, THUS
C          GENERATING AN EXTERNAL VIRTUAL SPACE.  WE SAVE
C          THE QUAMBO'S AND THIS RATHER ARBITRARY EXTERNAL SET.
C
C          -S- MATRIX AT WRK2, -Q- MATRIX AT WRK3
      CALL DAREAD(IDAF,IODA,WRK2,L2,12,0)
      CALL DAREAD(IDAF,IODA,WRK3,L3,45,0)
C           WE SEEK ORBITALS ORTHOGONAL TO THE -L1MBS- QUAMBO'S
      CALL ORTHO(WRK3,WRK2,WRK1,SCR,L1MBS,L0,L1,L2,L1)
      CALL TFSQB(WRK1,WRK3,SCR,L0,L1,L1)
      CALL DAWRIT(IDAF,IODA,WRK1,L3,470,0)
C
C       ----- TWO STEP GENERATION OF "SEMI-CANONICAL" ORBITALS -----
C          WE WANT TO DIAGONALIZE THE DIAGONAL BLOCKS OF F, ONCE
C          IN THE QUAMBO SPACE (OCCUPIED + VALENCE VIRTUAL) AND
C          ONCE IN THE EXTERNAL VIRTUAL SPACE.  WE SEEK TO HAVE
C          THE EIGENVALUES ORDERED IN EACH SPACE, BUT NOT TO MIX
C          VALENCE VIRTUALS AND EXTERNAL VIRTUALS.  WE ALSO WANT
C          THE DIAGONALIZATION TO EXPLOIT FULL POINT GROUP SYMMETRY.
C
C      --- CANONICALIZE BY FREE DIAGONALIZATION WITHIN EACH SPACE ---
C          TRANSFORM FOCK MATRIX TO QUAMBO+EXTERNAL ORBITAL SPACE,
C          THEN DIAGONALIZE THE SEPARATE BLOCKS, AND FINALLY
C          TRANSFORM THE VECTORS BACK TO THE ATOMIC BASIS SET.
C          VEC TO RECEIVE SEMI-CANONICAL QUAMBO/EXTERNAL SPACE
C
      CALL TFTRI(WRK3,FAO,WRK1,SCR,L0,L1,L1)
      CALL GLDIAG(L1,L1MBS,L1MBS,WRK3,SCR,EIG,VEC,IERR,IWRK)
      CALL TFSQB(VEC,WRK1,SCR,L1MBS,L1,L1)
      IJ=0
      DO I=L1MBS+1,L0
         IROW = (I*I-I)/2
         DO J=L1MBS+1,I
            IJ=IJ+1
            KL=IROW+J
            WRK3(IJ) = WRK3(KL)
         ENDDO
      ENDDO
      CALL GLDIAG(L1,L0-L1MBS,L0-L1MBS,WRK3,SCR,EIG(L1MBS+1),
     *            VEC(1,L1MBS+1),IERR,IWRK)
      CALL TFSQB(VEC(1,L1MBS+1),WRK1(L1MBS*L1+1),SCR,L0-L1MBS,L1,L1)
C
C      --- FULL SYMMETRY SEMI-CANONICALIZATION ---
C          THE ABOVE HAS PRODUCED ORBITALS ORDERED BY THEIR EIGENVALUE
C          WITHIN EACH SPACE, BUT LACKS CONTROL OF THE ORBITAL SYMMETRY.
C     A) BRING F TO PREVIOUSLY CANONICALIZED QUAMBO/EXTERNAL SPACE
C     B) ZAP ALL ELEMENTS CONNECTING THESE SPACES
C     C) TRANSFORM THIS MODIFIED FOCK OPERATOR BACK TO AO SPACE
C     D) DIAGONALIZE MODIFIED F WITH FULL Q-MATRIX SYMMETRY CONTROL
C     E) MAKE SURE THE SPACES DON'T GET MIXED UP BY OVERLAP MAXIMIZATION
C
C     THE RESULT OF THIS SHOULD BE THE FOLLOWING ORDER OF ORBITALS:
C            ORIGINAL OCCUPIED CANONICAL ORBITALS
C            VALENCE VIRTUAL ORBITALS
C            EXTERNAL VIRTUAL ORBITALS
C
      IF(NT.EQ.1) GO TO 500
C
C         -F IN SEMI-CANONICAL QUAMBO/EXTERNAL SPACE- AT WRK3
C         ONLY THE VIRTUAL VALENCE/EXTERNAL VIRTUAL BLOCK SHOULD
C         BE NONZERO, BUT ZAP OCCUPIED/EXTERNAL BLOCK TOO.
      CALL DAREAD(IDAF,IODA,FAO,L2,14,0)
      CALL TFTRI(WRK3,FAO,VEC,SCR,L0,L1,L1)
      DO I=L1MBS+1,L0
         I0 = (I*I-I)/2
         DO J=1,L1MBS
            WRK3(I0+J) = ZERO
         ENDDO
      ENDDO
C         -S- AT WRK2, -BLOCKED F IN MO SPACE- AT WRK3, SCRATCH AT WRK1
      CALL TFTRIB(FAO,WRK3,WRK2,VEC,WRK1,SCR,L0,L1,L2,L3)
C         STORE COPY OF ORBITALS AT -WRK1- FOR OVERLAP SELECTION BELOW
      CALL DCOPY(L0*L1,VEC,1,WRK1,1)
C         -Q- AT WRK3 USED TO SYMMETRY BLOCK DIAG STEP OF ZAPPED -FAO-
C         SCRATCH AT WRK2
      CALL DAREAD(IDAF,IODA,WRK3,L3,45,0)
      CALL TFTRI(WRK2,FAO,WRK3,SCR,L0,L1,L1)
      CALL SYMDIA(WRK2,VEC,EIG,SCR,IWRK,L0,L2,L1)
      CALL TFSQB(VEC,WRK3,SCR,L0,L1,L1)
C         ALIGN BY MAXIMUM OVERLAP TO FIRST STEP'S ORBITAL SPACE
C         -WRK1- IS ORIGINAL ORDER, -S- AT WRK2, SCRATCH AT WRK3
      CALL DAREAD(IDAF,IODA,WRK2,L2,12,0)
      CALL OVLSEL(VEC,EIG,WRK1,WRK2,WRK3,SCR,IWRK,SCR,L0,L1,L2)
C
C         ALL DONE, SAVE ORBITALS AND SEMI-CANONICAL EIGENVALUES.
C         NOTE THAT WE DO NOT NEED TO ASSIGN ORBITAL SYMMETRY LABELS,
C         WHICH IS TAKEN CARE OF IN THE HARTREE-FOCK DRIVER.
C
  500 CONTINUE
      CALL RETFM(NEED)
      IF(MASWRK) WRITE(IW,9090) NOCC,NVIR,1,NOCC,
     *           NOCC+1,NOCC+NVIRVAL,L1MBS,L0
      CALL DAWRIT(IDAF,IODA,VEC,L3,15,0)
      CALL DAWRIT(IDAF,IODA,EIG,L1,17,0)
      RETURN
C
 8000 FORMAT(/'FULLY ORTHONORMAL QUAMBO''S IN THE AO BASIS SET')
 8010 FORMAT(A5)
C
 9000 FORMAT(1X,'ABANDONING -VVO- GENERATION...'/
     *       1X,'THE NECESSARY ATOMIC AO SETS ONLY EXIST FOR H-AR,'/
     *       1X,'AND/OR PSEUDOPOTENTIALS MAY NOT BE USED.')
 9010 FORMAT(/1X,'FORMING VIRTUAL VALENCE ORBITALS...'/
     *       1X,'THE ATOMS POSSESS',I5,' CHEMICAL CORES, AND',I5,
     *          ' VALENCE ORBITALS.'/
     *       1X,'THIS MOLECULE HAS',I5,' OCCUPIED VALENCE ORBITALS,'/
     *       1X,'SO WE ARE SEEKING',I5,'  VIRTUAL VALENCE ORBITALS',
     *          ' FROM'/
     *       1X,'WITHIN THE MOLECULE''S VIRTUAL SPACE OF DIMENSION',I5)
 9020 FORMAT(1X,'TOO MANY SHELLS IN VVO ATOMIC MINIMAL BASIS, >',I8)
 9030 FORMAT(1X,'TOO MANY GAUSSIANS IN VVO ATOMIC MINIMAL BASIS, >',I8)
 9040 FORMAT(/1X, 'THERE SHOULD BE A CLEAR DROPOFF IN THE EIGENVALUES',
     *           ' OF -B- BETWEEN THE'/
     *           ' VALENCE VIRTUAL AND EXTERNAL VIRTUAL ORBITALS,',
     *           ' PLEASE CHECK THIS.')
 9050 FORMAT(1X,'LAST FEW EIGENVALUES KEPT AS VALENCE:')
 9060 FORMAT(1X,'FIRST FEW EIGENVALUES REJECTED AS EXTERNAL:')
 9070 FORMAT(1X,5F15.10)
 9080 FORMAT(/1X,'THE QUASIATOMIC MINIMAL-BASIS-SET ORBITALS',
     *          ' (QUAMBO''S) HAVE BEEN PUNCHED.'/
     *       1X,'THE',I5,' QUAMBO''S ARE MOLECULAR ORBITALS WHICH'/
     *       1X,'   A) ARE MUTUALLY ORTHONORMAL,'/
     *       1X,'   B) SPAN THE',I5,' DIMENSIONAL OCCUPIED ORBITAL',
     *          ' SPACE OF THIS MOLECULE'/
     *       1X,'      EXACTLY, AND',I5,
     *          ' ADDITIONAL VALENCE VIRTUAL ORBITALS,'/
     *       1X,'   C) ARE AS CLOSE AS POSSIBLE TO THE FREE ATOM AO''S',
     *          ' (IN THE SENSE OF '/
     *       1X,'      MAXIMUM OVERLAP), WHILE SATISFYING "A" AND "B".')
 9090 FORMAT(/1X,'VIRTUAL VALENCE ORBITAL GENERATION IS FINISHED.'/
     *       1X,'THE VIRTUAL VALENCE ORBITALS ARE A QUANTIFICATION OF',
     *          ' THE FRONTIER ORBITAL'/
     *       1X,'SPACE (I.E. LUMO, LUMO+1, ...) WITHIN THE HARTREE-',
     *          'FOCK APPROXIMATION.'/
     *       1X,'THE',I5,' OCCUPIED AND',I5,' VIRTUAL ORBITALS',
     *          ' PRINTED BELOW ARE'/
     *       1X,'     FIRST, THE OCCUPIED CANONICAL ORBITALS (',
     *           I5,'-',I5,')'/
     *       1X,'    SECOND, THE VIRTUAL VALENCE ORBITALS    (',
     *           I5,'-',I5,')'/
     *       1X,'     THIRD, THE EXTERNAL VIRTUAL ORBITALS   (',
     *           I5,'-',I5,')'/
     *       1X,'NOTE THAT THE VIRTUAL SPACE IS ONLY SEMI-',
     *          'CANONICALIZED, MEANING THAT'/
     *       1X,'THE FOCK OPERATOR IS DIAGONAL ONLY WITHIN ITS',
     *          ' VALENCE AND EXTERNAL'/
     *       1X,'VIRTUAL SUBBLOCKS.  THE VIRTUAL "EIGENVALUES" BELOW',
     *          ' ARE THESE DIAGONALS,'/
     *       1X,'WHICH ARE EXPECTATION VALUES OF THE FOCK OPERATOR,',
     *          ' NOT ITS EIGENVALUES.')
      END
C
C*MODULE VVOS    *DECK ETGTO
      SUBROUTINE ETGTO(NUCZ,DZ,CSINP,CPINP,IERR1,IERR2,LOC,
     *                 INTYP,NANGM,NBFS,MINF,MAXF,NSHLS,
     *                 EX,CS,CP,KSTART,KATOM,KTYPE,
     *                 KNG,KLOC,KMIN,KMAX,NGAUSS,NSHELL,MXGTOT,MXSH)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DONE,DZ,GOPARR,DSKWRK,MASWRK
C
      DIMENSION CSINP(MXGTOT),CPINP(MXGTOT),
     *          INTYP(*),NANGM(4),NBFS(4),MINF(4),MAXF(4),NSHLS(*),
     *          EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),
     *          KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *          KLOC(MXSH),KMIN(MXSH),KMAX(MXSH)
C         NEXT DIMENSIONS ARE FROM 18S,12P GTOS EXPANDING UP TO
C         FIVE SHELLS (1S,2S,2P,3S,3P) AND MAYBE FLOATING OUTER S,P
      DIMENSION EEX(30),COEF(80)
      DIMENSION IELIST(7),ICLIST(7),ITLIST(7),IGLIST(7)
C         4 HERE AND ABOVE IS BECAUSE PERIODIC TABLE HAS S,P,D,F OCCUP.
      CHARACTER*1 TYPE(4)
C
      PARAMETER (MXATM=2000)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA ZERO,ONE,PT5/0.0D+00,1.0D+00,0.5D+00/
      DATA PI32,TM6,TM10/5.56832799683170D+00,1.0D-06,1.0D-10/
      DATA TYPE/'S','P','D','F'/
C
C            ----- GET AN EVEN-TEMPERED BASIS SET -----
C     M.W.SCHMIDT, K.RUEDENBERG  J.CHEM.PHYS. 71, 3951-3962(1979)
C
C     THIS ROUTINE RETURNS EXPONENTS AND CONTRACTION COEFFICIENTS:
C
C     GTO EXPONENTS ARE THE OPTIMAL ALPHA,BETA VALUES FOR H-NE,
C     AND ARE FROM REGULARIZED ALPHA,BETA FOR NA-AR.  SIZES ARE
C     H-HE:8S  LI-BE:14S  B-NE:14S,7P  NA-MG:18S,9P  AL-AR:18S,12P
C     THESE ARE CHOSEN TO BE ABOUT ONE MILLIHARTREE OR SO AWAY FROM
C     THE EXACT ATOMIC SCF LIMITS (REPORTED IN THAT PAPER ALSO).
C
C     GENERAL CONTRACTION TO ATOMIC MINIMAL BASIS SET (DZ=.FALSE.),
C     OR FLOATING THE MOST DIFFUSE GTO IN EACH SYMMETRY (DZ=.TRUE.)
C
C     SPHERICAL DENSITY ATOMIC CALCULATIONS WILL GIVE ENERGIES
C         H -  FINAL ROHF ENERGY IS      -0.4999742933
C        HE -  FINAL RHF ENERGY IS       -2.8615070990
C        LI -  FINAL ROHF ENERGY IS      -7.4326894443
C        BE -  FINAL RHF ENERGY IS      -14.5729514471
C         B -  FINAL GVB ENERGY IS      -24.5289020592
C         C -  FINAL GVB ENERGY IS      -37.6883034560
C         N -  FINAL ROHF ENERGY IS     -54.4003685090
C         O -  FINAL GVB ENERGY IS      -74.8084211734
C         F -  FINAL GVB ENERGY IS      -99.4077883999
C        NE -  FINAL RHF ENERGY IS     -128.5447548872
C        NA -  FINAL ROHF ENERGY IS    -161.8583450857
C        MG -  FINAL RHF ENERGY IS     -199.6139973762
C        AL -  FINAL GVB ENERGY IS     -241.8758111063
C        SI -  FINAL GVB ENERGY IS     -288.8536108276
C         P -  FINAL ROHF ENERGY IS    -340.7179325784
C         S -  FINAL GVB ENERGY IS     -397.5038627824
C        CL -  FINAL GVB ENERGY IS     -459.4808442239
C        AR -  FINAL RHF ENERGY IS     -526.8160832671
C
      DO I=1,30
          EEX(I) = ZERO
      ENDDO
      DO I=1,80
         COEF(I) = ZERO
      ENDDO
C
C     ----- HYDROGEN TO HELIUM -----
C
      IF (NUCZ .GT. 2) GO TO 120
      CALL ETONE(EEX,COEF,NUCZ)
      GO TO 200
C
C     ----- LITHIUM TO NEON -----
C
  120 IF (NUCZ .GT. 10) GO TO 130
      CALL ETTWO(EEX,COEF,NUCZ)
      GO TO 200
C
C     ----- SODIUM TO ARGON -----
C
  130 IF (NUCZ .GT. 18) GO TO 190
      CALL ETTHREE(EEX,COEF,NUCZ)
      GO TO 200
C
C     ----- PAST ARGON DOES NOT EXIST -----
C
  190 CONTINUE
      CALL BERROR(3)
      STOP
C
C     ----- LOOP OVER EACH SHELL -----
C
  200 CONTINUE
      CALL ETSHELL(NUCZ,DZ,NPASS,IELIST,ICLIST,ITLIST,IGLIST)
C
      DONE  = .FALSE.
      IERR3 = 0
      IPASS = 0
  210 CONTINUE
      IPASS = IPASS+1
      IF(IPASS.GT.NPASS) DONE=.TRUE.
      IF(DONE. AND. IERR3.NE.0) CALL ABRT
      IF(DONE) RETURN
C
      NGE    = IELIST(IPASS)
      NGC    = ICLIST(IPASS)
      ITYP   = ITLIST(IPASS)
      IGAUSS = IGLIST(IPASS)
C
C     ----- DEFINE THE CURRENT SHELL -----
C
      NSHELL = NSHELL+1
      IF(NSHELL.GT.MXSH) THEN
         IERR1=1
         RETURN
      END IF
      NSHLS(NAT)     = NSHLS(NAT)+1
      KMIN(NSHELL)   = MINF(ITYP)
      KMAX(NSHELL)   = MAXF(ITYP)
      KSTART(NSHELL) = NGAUSS+1
      KATOM(NSHELL)  = NAT
      KTYPE(NSHELL)  = NANGM(ITYP)
      INTYP(NSHELL)  = ITYP
      KNG(NSHELL)    = IGAUSS
      KLOC(NSHELL)   = LOC+1
      NGAUSS         = NGAUSS+IGAUSS
      IF(NGAUSS.GT.MXGTOT) THEN
         IERR2=1
         RETURN
      END IF
      LOC = LOC+NBFS(ITYP)
      K1 = KSTART(NSHELL)
      K2 = K1+KNG(NSHELL)-1
      DO 440 I = 1,IGAUSS
         K = K1+I-1
         EX(K) = EEX(NGE+I)
         IF(ITYP.EQ.1) CSINP(K) = COEF(NGC+I)
         IF(ITYP.EQ.2) CPINP(K) = COEF(NGC+I)
         CS(K) = ZERO
         CP(K) = ZERO
         IF(ITYP.EQ.1) CS(K) = CSINP(K)
         IF(ITYP.EQ.2) CP(K) = CPINP(K)
  440 CONTINUE
C
C     ----- ALWAYS UNNORMALIZE PRIMITIVES -----
C
      DO 460 K = K1,K2
         EE = EX(K)+EX(K)
         FACS = PI32/(EE*SQRT(EE))
         FACP = PT5*FACS/EE
         IF(ITYP.EQ.1) CS(K) = CS(K)/SQRT(FACS)
         IF(ITYP.EQ.2) CP(K) = CP(K)/SQRT(FACP)
  460 CONTINUE
C
C     ----- IF(NORMF.EQ.0) NORMALIZE BASIS FUNCTIONS -----
C
      IF (NORMF .EQ. 1) GO TO 210
      FACS = ZERO
      FACP = ZERO
      DO 510 IG = K1,K2
         DO 500 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)
            IF (IG .EQ. JG) GO TO 480
               DUMS = DUMS+DUMS
               DUMP = DUMP+DUMP
  480       CONTINUE
            FACS = FACS+DUMS
            FACP = FACP+DUMP
  500    CONTINUE
  510 CONTINUE
C
      FAC=ZERO
      IF(ITYP.EQ.1 .AND. FACS.GT.TM10) FAC=ONE/SQRT(FACS*PI32)
      IF(ITYP.EQ.2 .AND. FACP.GT.TM10) FAC=ONE/SQRT(FACP*PI32)
C                                 VERIFY NORMALIZATION
      TNORM = ABS(ONE-FAC)
      IF(TNORM.LT.TM6) GO TO 520
C
      IF (MASWRK) THEN
         WRITE(IW,9000) FAC,NUCZ,IPASS,TYPE(ITYP)
         WRITE(IW,9010)  (EEX(NGE+K),K=1,IGAUSS)
         WRITE(IW,9020) (COEF(NGC+K),K=1,IGAUSS)
      ENDIF
      IERR3=IERR3+1
C
  520 CONTINUE
      DO 550 IG = K1,K2
         IF(ITYP.EQ.1) CS(IG) = FAC*CS(IG)
         IF(ITYP.EQ.2) CP(IG) = FAC*CP(IG)
         IF(ITYP.EQ.1) CSINP(IG) = FAC*CSINP(IG)
         IF(ITYP.EQ.2) CPINP(IG) = FAC*CPINP(IG)
  550 CONTINUE
      GO TO 210
C
 9000 FORMAT(/1X,'ERROR!!! NORMALIZATION FACTOR=',E16.8/
     *     1X,'FOR ATOM Z=',I3,' SHELL NO.',I3,' TYPE=',A1/
     *     1X,'CHECK BUILT IN ET-GTO EXPONENTS AND CONT. COEFS')
 9010 FORMAT(5F15.5)
 9020 FORMAT(5F15.8)
      END
C*MODULE VVOS    *DECK ETSHELL
      SUBROUTINE ETSHELL(NUCZ,DZ,NPASS,IELIST,ICLIST,ITLIST,IGLIST)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION IELIST(7),ICLIST(7),ITLIST(7),IGLIST(7)
      LOGICAL DZ
C
      CHARACTER*4 FLOAT
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C        FILL LISTS OF EXPONENT POINTER, COEFFICIENT POINTER,
C                      TYPE OF AO, AND NUMBER OF GAUSSIANS
C
      IF(NUCZ.GT. 2) GO TO 020
         NS =  8
         NP =  0
         NPASS=1
         ITLIST(1)=1
         IF(DZ) FLOAT='S   '
         GO TO 800
  020 IF(NUCZ.GT. 4) GO TO 030
         NS = 14
         NP =  0
         NPASS=2
         ITLIST(1)=1
         ITLIST(2)=1
         IF(DZ) FLOAT='S   '
         GO TO 800
  030 IF(NUCZ.GT.10) GO TO 040
         NS = 14
         NP =  7
         NPASS=3
         ITLIST(1)=1
         ITLIST(2)=1
         ITLIST(3)=2
         IF(DZ) FLOAT='SP  '
         GO TO 800
  040 IF(NUCZ.GT.12) GO TO 050
         NS = 18
         NP =  9
         NPASS=4
         ITLIST(1)=1
         ITLIST(2)=1
         ITLIST(3)=2
         ITLIST(4)=1
         IF(DZ) FLOAT='S   '
         GO TO 800
  050 IF(NUCZ.GT.18) GO TO 060
         NS = 18
         NP = 12
         NPASS=5
         ITLIST(1)=1
         ITLIST(2)=1
         ITLIST(3)=2
         ITLIST(4)=1
         ITLIST(5)=2
         IF(DZ) FLOAT='SP  '
         GO TO 800
C
C      THIS MISTAKE ALREADY TRAPPED FOR ELSEWHERE
C
  060 CONTINUE
      IF(MASWRK) WRITE(IW,*) 'IMPOSSIBLE MISTAKE IN ETSHELL'
      CALL ABRT
      STOP
C
  800 CONTINUE
      IELIST(1)=0
      ICLIST(1)=0
      IGLIST(1)=NS
      DO I=2,NPASS
         IF(ITLIST(I).EQ.1) THEN
            IELIST(I) = 0
            IGLIST(I) = NS
         END IF
         IF(ITLIST(I).EQ.2) THEN
            IELIST(I) = NS
            IGLIST(I) = NP
         END IF
         ICLIST(I) = ICLIST(I-1) + IGLIST(I-1)
      ENDDO
      IF(FLOAT.EQ.'S   ') THEN
         ITLIST(NPASS+1) = 1
         IGLIST(NPASS+1) = 1
         IELIST(NPASS+1) = 0
         ICLIST(NPASS+1) = ICLIST(NPASS) + IGLIST(NPASS)
         NPASS=NPASS+1
      END IF
      IF(FLOAT.EQ.'SP  ') THEN
         ITLIST(NPASS+1) = 1
         ITLIST(NPASS+2) = 2
         IGLIST(NPASS+1) = 1
         IGLIST(NPASS+2) = 1
         IELIST(NPASS+1) = 0
         IELIST(NPASS+2) = NS
         ICLIST(NPASS+1) = ICLIST(NPASS)   + IGLIST(NPASS)
         ICLIST(NPASS+2) = ICLIST(NPASS+1) + IGLIST(NPASS+1)
         NPASS=NPASS+2
      END IF
C
C-    WRITE(6,*) 'EXIT FROM ETSHELL ROUTINE, NPASS=',NPASS
C-    WRITE(6,899) 'IELIST',(IELIST(I),I=1,NPASS)
C-    WRITE(6,899) 'ICLIST',(ICLIST(I),I=1,NPASS)
C-    WRITE(6,899) 'ITLIST',(ITLIST(I),I=1,NPASS)
C-    WRITE(6,899) 'IGLIST',(IGLIST(I),I=1,NPASS)
C-899 FORMAT(A,'=',5I5)
C
      RETURN
      END
C*MODULE VVOS    *DECK ETONE
      SUBROUTINE ETONE(E,C,NUCZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION E(8),C(9)
C
      NS=8
      IF(NUCZ.EQ.2) GO TO 200
C
C           HYDROGEN
C
      ALPHA = 0.0304322D+00
      BETA  = 2.84371D+00
C         -1S-
      C( 1)= 2.78547941D-01
      C( 2)= 4.88609358D-01
      C( 3)= 2.54169291D-01
      C( 4)= 8.76801213D-02
      C( 5)= 2.50877638D-02
      C( 6)= 7.59416815D-03
      C( 7)= 1.49624961D-03
      C( 8)= 8.92506013D-04
      GO TO 5000
C
C           HELIUM
C
  200 CONTINUE
      ALPHA = 0.061317D+00
      BETA  = 3.06930D+00
C         -1S-
      C( 1)= 2.58538972D-01
      C( 2)= 4.68136547D-01
      C( 3)= 2.89598176D-01
      C( 4)= 1.12625711D-01
      C( 5)= 3.22554364D-02
      C( 6)= 8.81759572D-03
      C( 7)= 1.77431621D-03
      C( 8)= 8.16037658D-04
C
C        GENERATE THE EXPONENTS
C
 5000 CONTINUE
      E(1) = ALPHA*BETA
      DO I=2,NS
         E(I) = E(I-1) * BETA
      ENDDO
C
C        IN CASE WE ARE FLOATING THE OUTER GAUSSIAN
C
      C(9) = 1.0D+00
      RETURN
      END
C*MODULE VVOS    *DECK ETTWO
      SUBROUTINE ETTWO(E,C,NUCZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION E(21),C(37)
C
      IBR=NUCZ-2
      IF(IBR.LE.2) THEN
         NS=14
         NP=0
      ELSE
         NS=14
         NP=7
      ENDIF
      GO TO (100,200,300,400,500,600,700,800), IBR
C
C           LITHIUM
C
  100 CONTINUE
      ALPHAS = 0.010167D+00
      BETAS  = 2.57870D+00
      ALPHAP = 0.0D+00
      BETAP  = 0.0D+00
C         -1S-
      C( 1)=-7.29556274D-04
      C( 2)= 3.72933364D-03
      C( 3)=-1.45835583D-03
      C( 4)= 1.74570321D-01
      C( 5)= 4.12462766D-01
      C( 6)= 3.24701318D-01
      C( 7)= 1.59239563D-01
      C( 8)= 6.09839168D-02
      C( 9)= 2.02383744D-02
      C(10)= 6.50231488D-03
      C(11)= 1.91557863D-03
      C(12)= 6.87713606D-04
      C(13)= 1.16324310D-04
      C(14)= 1.10197460D-04
C         -2S-
      C(15)= 4.56906084D-01
      C(16)= 6.08640023D-01
      C(17)= 4.77727591D-02
      C(18)=-1.21468704D-01
      C(19)=-1.02607549D-01
      C(20)=-5.85732292D-02
      C(21)=-2.66406189D-02
      C(22)=-9.74138553D-03
      C(23)=-3.21171832D-03
      C(24)=-1.01465541D-03
      C(25)=-3.01728226D-04
      C(26)=-1.06865730D-04
      C(27)=-1.83748435D-05
      C(28)=-1.71804556D-05
      GO TO 5000
C
C           BERYLLIUM
C
  200 CONTINUE
      ALPHAS = 0.017502D+00
      BETAS  = 2.58803D+00
      ALPHAP = 0.0D+00
      BETAP  = 0.0D+00
C         -1S-
      C( 1)=-8.32062811D-04
      C( 2)= 2.87338821D-03
      C( 3)=-3.90449716D-03
      C( 4)= 1.40044621D-01
      C( 5)= 4.12304599D-01
      C( 6)= 3.45796460D-01
      C( 7)= 1.67498022D-01
      C( 8)= 6.37016457D-02
      C( 9)= 2.08846484D-02
      C(10)= 6.70317405D-03
      C(11)= 1.95650359D-03
      C(12)= 7.00733918D-04
      C(13)= 1.18699193D-04
      C(14)= 1.10772010D-04
C         -2S-
      C(15)= 3.43871894D-01
      C(16)= 5.98134289D-01
      C(17)= 2.07436470D-01
      C(18)=-1.05687931D-01
      C(19)=-1.34940265D-01
      C(20)=-7.55299312D-02
      C(21)=-3.32753335D-02
      C(22)=-1.18934244D-02
      C(23)=-3.88866372D-03
      C(24)=-1.21585860D-03
      C(25)=-3.61380098D-04
      C(26)=-1.26412191D-04
      C(27)=-2.20513869D-05
      C(28)=-2.01071534D-05
      GO TO 5000
C
C           BORON
C
  300 CONTINUE
      ALPHAS = 0.029291D+00
      BETAS  = 2.58041D+00
      ALPHAP = 0.021928D+00
      BETAP  = 2.80107D+00
C         -1S-
      C( 1)=-7.51750475D-04
      C( 2)= 3.09099384D-03
      C( 3)=-2.78589641D-03
      C( 4)= 1.41804573D-01
      C( 5)= 4.15219672D-01
      C( 6)= 3.44919426D-01
      C( 7)= 1.63805046D-01
      C( 8)= 6.16166843D-02
      C( 9)= 2.01312598D-02
      C(10)= 6.47919925D-03
      C(11)= 1.89252354D-03
      C(12)= 6.83036765D-04
      C(13)= 1.14603380D-04
      C(14)= 1.08957043D-04
C         -2S-
      C(15)= 3.05077843D-01
      C(16)= 5.90329023D-01
      C(17)= 2.61905241D-01
      C(18)=-9.89170234D-02
      C(19)=-1.51949624D-01
      C(20)=-8.47272362D-02
      C(21)=-3.62302474D-02
      C(22)=-1.27734987D-02
      C(23)=-4.15524044D-03
      C(24)=-1.30279436D-03
      C(25)=-3.87264329D-04
      C(26)=-1.36568885D-04
      C(27)=-2.35838704D-05
      C(28)=-2.19163253D-05
C         -2P-
      C(29)= 2.36473469D-01
      C(30)= 4.73461866D-01
      C(31)= 3.39857244D-01
      C(32)= 1.34359035D-01
      C(33)= 3.80827585D-02
      C(34)= 8.86824791D-03
      C(35)= 2.49613549D-03
      GO TO 5000
C
C           CARBON
C
  400 CONTINUE
      ALPHAS = 0.043391D+00
      BETAS  = 2.57899D+00
      ALPHAP = 0.032290D+00
      BETAP  = 2.81511D+00
C         -1S-
      C( 1)=-6.77891849D-04
      C( 2)= 3.16694226D-03
      C( 3)=-1.91224170D-03
      C( 4)= 1.40383478D-01
      C( 5)= 4.16044863D-01
      C( 6)= 3.46268829D-01
      C( 7)= 1.62633679D-01
      C( 8)= 6.07386990D-02
      C( 9)= 1.97742578D-02
      C(10)= 6.36639261D-03
      C(11)= 1.85791680D-03
      C(12)= 6.71969050D-04
      C(13)= 1.12417139D-04
      C(14)= 1.07343121D-04
C         -2S-
      C(15)= 2.86217712D-01
      C(16)= 5.79928125D-01
      C(17)= 2.92836401D-01
      C(18)=-8.98781103D-02
      C(19)=-1.61868757D-01
      C(20)=-9.07202710D-02
      C(21)=-3.80900959D-02
      C(22)=-1.33060046D-02
      C(23)=-4.30857739D-03
      C(24)=-1.35123146D-03
      C(25)=-4.01167127D-04
      C(26)=-1.41800712D-04
      C(27)=-2.44083143D-05
      C(28)=-2.27858047D-05
C         -2P-
      C(29)= 2.06607795D-01
      C(30)= 4.55793793D-01
      C(31)= 3.65759456D-01
      C(32)= 1.56210647D-01
      C(33)= 4.46494975D-02
      C(34)= 9.97314134D-03
      C(35)= 2.69592178D-03
      GO TO 5000
C
C           NITROGEN
C
  500 CONTINUE
      ALPHAS = 0.060829D+00
      BETAS  = 2.57309D+00
      ALPHAP = 0.046701D+00
      BETAP  = 2.81544D+00
C         -1S-
      C( 1)=-2.35953751D-03
      C( 2)=-3.10439370D-04
      C( 3)=-3.00749064D-03
      C( 4)= 1.40467696D-01
      C( 5)= 4.16709114D-01
      C( 6)= 3.47340225D-01
      C( 7)= 1.62298128D-01
      C( 8)= 6.04634815D-02
      C( 9)= 1.96908729D-02
      C(10)= 6.35762183D-03
      C(11)= 1.85766295D-03
      C(12)= 6.75771050D-04
      C(13)= 1.12251348D-04
      C(14)= 1.08728458D-04
C         -2S-
      C(15)= 2.81652533D-01
      C(16)= 5.68964214D-01
      C(17)= 3.08574614D-01
      C(18)=-8.24456453D-02
      C(19)=-1.65103241D-01
      C(20)=-9.25844872D-02
      C(21)=-3.83425399D-02
      C(22)=-1.33290246D-02
      C(23)=-4.30911285D-03
      C(24)=-1.35559342D-03
      C(25)=-4.02550112D-04
      C(26)=-1.43231592D-04
      C(27)=-2.44521149D-05
      C(28)=-2.31743715D-05
C         -2P-
      C(29)= 2.06671781D-01
      C(30)= 4.43627998D-01
      C(31)= 3.71713817D-01
      C(32)= 1.64190671D-01
      C(33)= 4.71480615D-02
      C(34)= 1.02400330D-02
      C(35)= 2.70620542D-03
      GO TO 5000
C
C           OXYGEN
C
  600 CONTINUE
      ALPHAS = 0.080317D+00
      BETAS  = 2.57088D+00
      ALPHAP = 0.054772D+00
      BETAP  = 2.85531D+00
C         -1S-
      C( 1)=-5.46447458D-04
      C( 2)= 3.16835882D-03
      C( 3)=-6.93321355D-04
      C( 4)= 1.37402996D-01
      C( 5)= 4.14838002D-01
      C( 6)= 3.48691078D-01
      C( 7)= 1.62580253D-01
      C( 8)= 6.04707128D-02
      C( 9)= 1.96824379D-02
      C(10)= 6.36296563D-03
      C(11)= 1.85932316D-03
      C(12)= 6.78107985D-04
      C(13)= 1.12273370D-04
      C(14)= 1.09384895D-04
C         -2S-
      C(15)= 2.67229609D-01
      C(16)= 5.63077851D-01
      C(17)= 3.28757979D-01
      C(18)=-7.64700218D-02
      C(19)=-1.73325650D-01
      C(20)=-9.88222677D-02
      C(21)=-4.07278039D-02
      C(22)=-1.41474229D-02
      C(23)=-4.56540797D-03
      C(24)=-1.44002060D-03
      C(25)=-4.26973251D-04
      C(26)=-1.52565143D-04
      C(27)=-2.59039126D-05
      C(28)=-2.47353092D-05
C         -2P-
      C(29)= 1.95485469D-01
      C(30)= 4.21336327D-01
      C(31)= 3.86683918D-01
      C(32)= 1.87487017D-01
      C(33)= 5.47811310D-02
      C(34)= 1.15836175D-02
      C(35)= 2.93577527D-03
      GO TO 5000
C
C           FLUORINE
C
  700 CONTINUE
      ALPHAS = 0.102117D+00
      BETAS  = 2.56967D+00
      ALPHAP = 0.065739D+00
      BETAP  = 2.87614D+00
C         -1S-
      C( 1)=-4.87577685D-04
      C( 2)= 3.11283027D-03
      C( 3)=-3.45388158D-04
      C( 4)= 1.34668231D-01
      C( 5)= 4.13750035D-01
      C( 6)= 3.50764600D-01
      C( 7)= 1.63399997D-01
      C( 8)= 6.07070524D-02
      C( 9)= 1.97490310D-02
      C(10)= 6.38938129D-03
      C(11)= 1.86684354D-03
      C(12)= 6.81895494D-04
      C(13)= 1.12677736D-04
      C(14)= 1.10145916D-04
C         -2S-
      C(15)= 2.57310911D-01
      C(16)= 5.57447726D-01
      C(17)= 3.43251149D-01
      C(18)=-6.98130905D-02
      C(19)=-1.77303422D-01
      C(20)=-1.02068035D-01
      C(21)=-4.18409023D-02
      C(22)=-1.45058338D-02
      C(23)=-4.67322314D-03
      C(24)=-1.47589369D-03
      C(25)=-4.37182995D-04
      C(26)=-1.56579848D-04
      C(27)=-2.65034847D-05
      C(28)=-2.54138755D-05
C         -2P-
      C(29)= 1.86334154D-01
      C(30)= 4.07371202D-01
      C(31)= 3.95329733D-01
      C(32)= 2.02798538D-01
      C(33)= 6.03786615D-02
      C(34)= 1.25455726D-02
      C(35)= 3.10771124D-03
      GO TO 5000
C
C           NEON
C
  800 CONTINUE
      ALPHAS = 0.127936D+00
      BETAS  = 2.56563D+00
      ALPHAP = 0.079539D+00
      BETAP  = 2.89200D+00
C         -1S-
      C( 1)=-4.50483068D-04
      C( 2)= 3.11348476D-03
      C( 3)= 4.36818623D-05
      C( 4)= 1.34299535D-01
      C( 5)= 4.13072882D-01
      C( 6)= 3.51120594D-01
      C( 7)= 1.63439881D-01
      C( 8)= 6.07117712D-02
      C( 9)= 1.97712069D-02
      C(10)= 6.40967134D-03
      C(11)= 1.87484231D-03
      C(12)= 6.87458175D-04
      C(13)= 1.13052306D-04
      C(14)= 1.11600678D-04
C         -2S-
      C(15)= 2.55997387D-01
      C(16)= 5.51817518D-01
      C(17)= 3.49753139D-01
      C(18)=-6.62818503D-02
      C(19)=-1.79682749D-01
      C(20)=-1.04118761D-01
      C(21)=-4.25120285D-02
      C(22)=-1.47349754D-02
      C(23)=-4.74592434D-03
      C(24)=-1.50324875D-03
      C(25)=-4.45227280D-04
      C(26)=-1.60272951D-04
      C(27)=-2.69542751D-05
      C(28)=-2.61332647D-05
C         -2P-
      C(29)= 1.81887066D-01
      C(30)= 3.99529131D-01
      C(31)= 4.00021287D-01
      C(32)= 2.11238040D-01
      C(33)= 6.33759975D-02
      C(34)= 1.29660837D-02
      C(35)= 3.14458087D-03
C
C        GENERATE THE EXPONENTS
C
 5000 CONTINUE
      E(1) = ALPHAS*BETAS
      DO I=2,NS
         E(I) = E(I-1) * BETAS
      ENDDO
      E(NS+1) = ALPHAP*BETAP
      DO I=2,NP
         E(I+NS) = E(I+NS-1) * BETAP
      ENDDO
C
C        IN CASE WE ARE FLOATING THE OUTER GAUSSIANS
C
      IF(NUCZ.LE.4) THEN
         C(29) = 1.0D+00
      ELSE
         C(36) = 1.0D+00
         C(37) = 1.0D+00
      END IF
      RETURN
      END
C*MODULE VVOS    *DECK ETTHREE
      SUBROUTINE ETTHREE(E,C,NUCZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION E(30),C(80)
C
      IBR=NUCZ-10
      IF(IBR.LE.2) THEN
         NS=18
         NP= 9
      ELSE
         NS=18
         NP=12
      ENDIF
      GO TO (100,200,300,400,500,600,700,800), IBR
C
C           SODIUM
C
  100 CONTINUE
      ALPHAS = 9.829283D-03
      BETAS  = 2.490586D+00
      ALPHAP = 1.115268D-01
      BETAP  = 2.526571D+00
C         -1S-
      C( 1)= 1.84295524D-05
      C( 2)=-6.91485856D-05
      C( 3)= 1.58616496D-04
      C( 4)=-2.44260865D-04
      C( 5)= 1.52174325D-03
      C( 6)= 1.17336136D-03
      C( 7)= 8.89556694D-02
      C( 8)= 3.68360801D-01
      C( 9)= 3.78250276D-01
      C(10)= 1.96767810D-01
      C(11)= 7.85275924D-02
      C(12)= 2.72154075D-02
      C(13)= 9.03852939D-03
      C(14)= 2.95462201D-03
      C(15)= 9.00388808D-04
      C(16)= 3.46650172D-04
      C(17)= 5.29159543D-05
      C(18)= 6.22831358D-05
C         -2S-
      C(19)=-7.76271392D-04
      C(20)= 3.51969796D-03
      C(21)=-5.65958968D-03
      C(22)= 1.67055626D-01
      C(23)= 5.50601638D-01
      C(24)= 4.16026074D-01
      C(25)=-1.14712906D-02
      C(26)=-1.84090557D-01
      C(27)=-1.23937468D-01
      C(28)=-5.45055029D-02
      C(29)=-2.01018427D-02
      C(30)=-6.81312923D-03
      C(31)=-2.21541789D-03
      C(32)=-7.27771992D-04
      C(33)=-2.18783526D-04
      C(34)=-8.51812585D-05
      C(35)=-1.27775065D-05
      C(36)=-1.52405315D-05
C         -2P-
      C(37)= 9.39335927D-02
      C(38)= 3.21096694D-01
      C(39)= 3.96698886D-01
      C(40)= 2.76010656D-01
      C(41)= 1.18682254D-01
      C(42)= 3.48793987D-02
      C(43)= 8.98823348D-03
      C(44)= 1.72837476D-03
      C(45)= 5.62829087D-04
C         -3S-
      C(46)= 4.99571681D-01
      C(47)= 5.93164405D-01
      C(48)= 2.58846077D-02
      C(49)=-1.45189986D-01
      C(50)=-1.30965138D-01
      C(51)=-7.37744668D-02
      C(52)= 1.32610037D-03
      C(53)= 2.89990233D-02
      C(54)= 1.89224257D-02
      C(55)= 8.26397024D-03
      C(56)= 3.02469697D-03
      C(57)= 1.02702617D-03
      C(58)= 3.32385038D-04
      C(59)= 1.09652556D-04
      C(60)= 3.27861986D-05
      C(61)= 1.28345294D-05
      C(62)= 1.91013312D-06
      C(63)= 2.29238790D-06
      GO TO 5000
C
C           MAGNESIUM
C
  200 CONTINUE
      ALPHAS = 1.453475D-02
      BETAS  = 2.464371D+00
      ALPHAP = 1.644377D-01
      BETAP  = 2.483089D+00
C         -1S-
      C( 1)= 9.51104928D-05
      C( 2)=-3.31543241D-04
      C( 3)= 7.11517704D-04
      C( 4)=-1.19834019D-03
      C( 5)= 3.10919122D-03
      C( 6)= 1.23106656D-03
      C( 7)= 1.24815898D-01
      C( 8)= 3.89532074D-01
      C( 9)= 3.53355621D-01
      C(10)= 1.76376786D-01
      C(11)= 6.95857347D-02
      C(12)= 2.42326563D-02
      C(13)= 8.11332108D-03
      C(14)= 2.69339021D-03
      C(15)= 8.24537433D-04
      C(16)= 3.27397642D-04
      C(17)= 4.74315261D-05
      C(18)= 6.07679547D-05
C         -2S-
      C(19)=-1.13783033D-03
      C(20)= 4.09545780D-03
      C(21)=-6.53269113D-03
      C(22)= 1.95198386D-01
      C(23)= 5.88200544D-01
      C(24)= 3.72677400D-01
      C(25)=-5.92541897D-02
      C(26)=-1.87989662D-01
      C(27)=-1.16363354D-01
      C(28)=-5.00495794D-02
      C(29)=-1.83533575D-02
      C(30)=-6.26948888D-03
      C(31)=-2.05586699D-03
      C(32)=-6.86547439D-04
      C(33)=-2.07192495D-04
      C(34)=-8.32615352D-05
      C(35)=-1.18348951D-05
      C(36)=-1.53860002D-05
C         -2P-
      C(37)= 9.18456923D-02
      C(38)= 3.30709339D-01
      C(39)= 4.01662845D-01
      C(40)= 2.63523357D-01
      C(41)= 1.09152615D-01
      C(42)= 3.17079546D-02
      C(43)= 8.26724283D-03
      C(44)= 1.58770689D-03
      C(45)= 5.48400537D-04
C         -3S-
      C(46)= 3.95634766D-01
      C(47)= 6.15732013D-01
      C(48)= 1.64776758D-01
      C(49)=-1.83049695D-01
      C(50)=-1.87828525D-01
      C(51)=-8.47240834D-02
      C(52)= 1.16685412D-02
      C(53)= 3.81270176D-02
      C(54)= 2.27537301D-02
      C(55)= 9.74562053D-03
      C(56)= 3.53508559D-03
      C(57)= 1.21359255D-03
      C(58)= 3.94611927D-04
      C(59)= 1.32883477D-04
      C(60)= 3.96927998D-05
      C(61)= 1.61173623D-05
      C(62)= 2.25584657D-06
      C(63)= 2.96866028D-06
      GO TO 5000
C
C           ALUMINUM
C
  300 CONTINUE
      ALPHAS = 2.392865D-02
      BETAS  = 2.434760D+00
      ALPHAP = 9.540427D-03
      BETAP  = 2.555805D+00
C         -1S-
      C( 1)= 1.82921603D-04
      C( 2)=-6.40215464D-04
      C( 3)= 1.33910138D-03
      C( 4)=-2.14806231D-03
      C( 5)= 4.70750056D-03
      C( 6)= 9.69030627D-03
      C( 7)= 2.01280358D-01
      C( 8)= 4.10008452D-01
      C( 9)= 3.03885978D-01
      C(10)= 1.41831537D-01
      C(11)= 5.44303203D-02
      C(12)= 1.90451750D-02
      C(13)= 6.39422633D-03
      C(14)= 2.16997761D-03
      C(15)= 6.63043970D-04
      C(16)= 2.74794024D-04
      C(17)= 3.67944377D-05
      C(18)= 5.28129236D-05
C         -2S-
      C(19)=-1.76010692D-03
      C(20)= 6.38708887D-03
      C(21)=-1.71297331D-03
      C(22)= 3.04506458D-01
      C(23)= 6.14107816D-01
      C(24)= 2.49271849D-01
      C(25)=-1.30535154D-01
      C(26)=-1.80454861D-01
      C(27)=-9.78551117D-02
      C(28)=-4.06016405D-02
      C(29)=-1.46492349D-02
      C(30)=-5.05380465D-03
      C(31)=-1.66444292D-03
      C(32)=-5.68851914D-04
      C(33)=-1.71314653D-04
      C(34)=-7.18910434D-05
      C(35)=-9.43185916D-06
      C(36)=-1.37553286D-05
C         -2P-
      C(37)= 4.19041372D-04
      C(38)=-9.05020149D-04
      C(39)= 3.01320536D-03
      C(40)= 3.31873496D-02
      C(41)= 2.78898528D-01
      C(42)= 4.35727941D-01
      C(43)= 3.05103398D-01
      C(44)= 1.27627802D-01
      C(45)= 3.60577208D-02
      C(46)= 8.89576826D-03
      C(47)= 1.66380228D-03
      C(48)= 5.18945801D-04
C         -3S-
      C(49)= 3.99091023D-01
      C(50)= 6.29832451D-01
      C(51)= 1.72625254D-01
      C(52)=-2.54601846D-01
      C(53)=-2.15976240D-01
      C(54)=-6.77813858D-02
      C(55)= 3.24868809D-02
      C(56)= 4.35258513D-02
      C(57)= 2.29412277D-02
      C(58)= 9.45610861D-03
      C(59)= 3.38756276D-03
      C(60)= 1.17051611D-03
      C(61)= 3.83798039D-04
      C(62)= 1.31682366D-04
      C(63)= 3.94574508D-05
      C(64)= 1.66407142D-05
      C(65)= 2.16632820D-06
      C(66)= 3.17893798D-06
C         -3P-
      C(67)= 1.12544997D-01
      C(68)= 4.83567560D-01
      C(69)= 4.37010409D-01
      C(70)= 1.25811071D-01
      C(71)=-6.06740437D-02
      C(72)=-8.54431097D-02
      C(73)=-5.76406082D-02
      C(74)=-2.36159708D-02
      C(75)=-6.55642601D-03
      C(76)=-1.60624413D-03
      C(77)=-2.99494373D-04
      C(78)=-9.33337552D-05
      GO TO 5000
C
C           SILICON
C
  400 CONTINUE
      ALPHAS = 3.336845D-02
      BETAS  = 2.418537D+00
      ALPHAP = 2.246645D-02
      BETAP  = 2.443551D+00
C         -1S-
      C( 1)= 1.76315111D-04
      C( 2)=-6.23296388D-04
      C( 3)= 1.27344907D-03
      C( 4)=-1.93916104D-03
      C( 5)= 4.34644658D-03
      C( 6)= 2.12004217D-02
      C( 7)= 2.43886724D-01
      C( 8)= 4.09222255D-01
      C( 9)= 2.76903707D-01
      C(10)= 1.25596359D-01
      C(11)= 4.75868402D-02
      C(12)= 1.67335018D-02
      C(13)= 5.62359529D-03
      C(14)= 1.93479034D-03
      C(15)= 5.89488197D-04
      C(16)= 2.50739183D-04
      C(17)= 3.18897731D-05
      C(18)= 4.90964745D-05
C         -2S-
      C(19)=-1.93307676D-03
      C(20)= 7.35442611D-03
      C(21)= 4.27840082D-03
      C(22)= 3.60970377D-01
      C(23)= 6.15159584D-01
      C(24)= 1.83233397D-01
      C(25)=-1.59519967D-01
      C(26)=-1.74270300D-01
      C(27)=-8.94103640D-02
      C(28)=-3.64638404D-02
      C(29)=-1.30791051D-02
      C(30)=-4.53596105D-03
      C(31)=-1.49972771D-03
      C(32)=-5.18591513D-04
      C(33)=-1.56175328D-04
      C(34)=-6.70760241D-05
      C(35)=-8.39159975D-06
      C(36)=-1.30874527D-05
C         -2P-
      C(37)=-3.21476543D-04
      C(38)= 1.88559789D-03
      C(39)= 2.13109218D-03
      C(40)= 1.11297472D-01
      C(41)= 3.73669790D-01
      C(42)= 3.95073439D-01
      C(43)= 2.24933142D-01
      C(44)= 8.34050497D-02
      C(45)= 2.31258208D-02
      C(46)= 5.91084670D-03
      C(47)= 1.12355276D-03
      C(48)= 4.05476948D-04
C         -3S-
      C(49)= 3.84365606D-01
      C(50)= 6.44187329D-01
      C(51)= 1.92141878D-01
      C(52)=-2.98690395D-01
      C(53)=-2.36445899D-01
      C(54)=-5.51012485D-02
      C(55)= 4.42109443D-02
      C(56)= 4.65419083D-02
      C(57)= 2.32481658D-02
      C(58)= 9.41757999D-03
      C(59)= 3.35578052D-03
      C(60)= 1.16515234D-03
      C(61)= 3.83736817D-04
      C(62)= 1.33132392D-04
      C(63)= 3.99188536D-05
      C(64)= 1.72178630D-05
      C(65)= 2.13941904D-06
      C(66)= 3.35489933D-06
C         -3P-
      C(67)= 2.55168574D-01
      C(68)= 5.00630509D-01
      C(69)= 3.49675355D-01
      C(70)= 3.73720506D-02
      C(71)=-9.62282726D-02
      C(72)=-9.04739191D-02
      C(73)=-5.07163363D-02
      C(74)=-1.82036236D-02
      C(75)=-5.01971708D-03
      C(76)=-1.25942127D-03
      C(77)=-2.43338699D-04
      C(78)=-8.61066678D-05
      GO TO 5000
C
C           PHOSPHORUS
C
  500 CONTINUE
      ALPHAS = 4.340876D-02
      BETAS  = 2.403573D+00
      ALPHAP = 3.076660D-02
      BETAP  = 2.410779D+00
C         -1S-
      C( 1)= 1.84330099D-04
      C( 2)=-4.77111036D-04
      C( 3)= 1.08154764D-03
      C( 4)=-1.62011519D-03
      C( 5)= 3.59260401D-03
      C( 6)= 3.09244746D-02
      C( 7)= 2.67459444D-01
      C( 8)= 4.04463903D-01
      C( 9)= 2.61621178D-01
      C(10)= 1.17276112D-01
      C(11)= 4.42885547D-02
      C(12)= 1.56739763D-02
      C(13)= 5.28637474D-03
      C(14)= 1.84015463D-03
      C(15)= 5.60152640D-04
      C(16)= 2.43763131D-04
      C(17)= 2.95229314D-05
      C(18)= 4.85869978D-05
C         -2S-
      C(19)=-2.98946064D-03
      C(20)= 6.27209506D-03
      C(21)= 7.97228212D-03
      C(22)= 3.92527441D-01
      C(23)= 6.14535292D-01
      C(24)= 1.45290016D-01
      C(25)=-1.74123438D-01
      C(26)=-1.70824460D-01
      C(27)=-8.53839452D-02
      C(28)=-3.45830653D-02
      C(29)=-1.24098968D-02
      C(30)=-4.32962145D-03
      C(31)=-1.43942102D-03
      C(32)=-5.02753014D-04
      C(33)=-1.51594135D-04
      C(34)=-6.64686757D-05
      C(35)=-7.94303749D-06
      C(36)=-1.32106755D-05
C         -2P-
      C(37)=-4.02211487D-04
      C(38)= 1.78886775D-03
      C(39)= 2.57826540D-03
      C(40)= 1.17964187D-01
      C(41)= 3.81614412D-01
      C(42)= 3.89880344D-01
      C(43)= 2.15172682D-01
      C(44)= 7.88617719D-02
      C(45)= 2.19388600D-02
      C(46)= 5.69938108D-03
      C(47)= 1.08717717D-03
      C(48)= 4.12414147D-04
C         -3S-
      C(49)= 3.68459686D-01
      C(50)= 6.53131508D-01
      C(51)= 2.15189889D-01
      C(52)=-3.26925800D-01
      C(53)=-2.52709754D-01
      C(54)=-4.60522195D-02
      C(55)= 5.12279304D-02
      C(56)= 4.84288498D-02
      C(57)= 2.35330608D-02
      C(58)= 9.47677720D-03
      C(59)= 3.37402174D-03
      C(60)= 1.17995007D-03
      C(61)= 3.90176858D-04
      C(62)= 1.36954321D-04
      C(63)= 4.10357062D-05
      C(64)= 1.81045523D-05
      C(65)= 2.14160500D-06
      C(66)= 3.59127671D-06
C         -3P-
      C(67)= 2.30769800D-01
      C(68)= 4.90530633D-01
      C(69)= 3.77388620D-01
      C(70)= 4.59711706D-02
      C(71)=-1.08884423D-01
      C(72)=-1.00087890D-01
      C(73)=-5.41873817D-02
      C(74)=-1.91659027D-02
      C(75)=-5.30378214D-03
      C(76)=-1.35018780D-03
      C(77)=-2.62292734D-04
      C(78)=-9.73778566D-05
      GO TO 5000
C
C           SULFUR
C
  600 CONTINUE
      ALPHAS = 5.481334D-02
      BETAS  = 2.388220D+00
      ALPHAP = 3.572023D-02
      BETAP  = 2.400897D+00
C         -1S-
      C( 1)= 1.19033643D-04
      C( 2)=-4.39215001D-04
      C( 3)= 8.30269610D-04
      C( 4)=-1.11714029D-03
      C( 5)= 3.08817528D-03
      C( 6)= 3.95486573D-02
      C( 7)= 2.83078881D-01
      C( 8)= 3.98919917D-01
      C( 9)= 2.51014812D-01
      C(10)= 1.11951756D-01
      C(11)= 4.23185873D-02
      C(12)= 1.50870580D-02
      C(13)= 5.11470734D-03
      C(14)= 1.80084567D-03
      C(15)= 5.48173559D-04
      C(16)= 2.44113953D-04
      C(17)= 2.80299594D-05
      C(18)= 4.95694890D-05
C         -2S-
      C(19)=-1.96130999D-03
      C(20)= 8.52889331D-03
      C(21)= 1.25393487D-02
      C(22)= 4.14035380D-01
      C(23)= 6.10924767D-01
      C(24)= 1.17318575D-01
      C(25)=-1.83333374D-01
      C(26)=-1.68258638D-01
      C(27)=-8.28534488D-02
      C(28)=-3.34891518D-02
      C(29)=-1.20567930D-02
      C(30)=-4.23571168D-03
      C(31)=-1.41720765D-03
      C(32)=-5.00113186D-04
      C(33)=-1.51005674D-04
      C(34)=-6.76579102D-05
      C(35)=-7.68042516D-06
      C(36)=-1.37053471D-05
C         -2P-
      C(37)=-3.91510033D-05
      C(38)= 1.44961795D-03
      C(39)= 2.33617957D-03
      C(40)= 9.82346305D-02
      C(41)= 3.70817374D-01
      C(42)= 4.01205338D-01
      C(43)= 2.24336645D-01
      C(44)= 8.33139090D-02
      C(45)= 2.32697094D-02
      C(46)= 6.10666792D-03
      C(47)= 1.16244802D-03
      C(48)= 4.49994519D-04
C         -3S-
      C(49)= 3.56611583D-01
      C(50)= 6.57653646D-01
      C(51)= 2.37169823D-01
      C(52)=-3.53350211D-01
      C(53)=-2.69330236D-01
      C(54)=-3.93723801D-02
      C(55)= 5.76414360D-02
      C(56)= 5.09832196D-02
      C(57)= 2.44003607D-02
      C(58)= 9.81164219D-03
      C(59)= 3.50271486D-03
      C(60)= 1.23425077D-03
      C(61)= 4.10445865D-04
      C(62)= 1.45674484D-04
      C(63)= 4.36668579D-05
      C(64)= 1.97047091D-05
      C(65)= 2.21006305D-06
      C(66)= 3.98265583D-06
C         -3P-
      C(67)= 2.05853411D-01
      C(68)= 4.52591378D-01
      C(69)= 4.13470723D-01
      C(70)= 9.22401174D-02
      C(71)=-1.12872114D-01
      C(72)=-1.11811324D-01
      C(73)=-6.08088989D-02
      C(74)=-2.17783609D-02
      C(75)=-6.02983392D-03
      C(76)=-1.55266967D-03
      C(77)=-2.99956773D-04
      C(78)=-1.13975996D-04
      GO TO 5000
C
C           CHLORINE
C
  700 CONTINUE
      ALPHAS = 6.679696D-02
      BETAS  = 2.380215D+00
      ALPHAP = 4.122936D-02
      BETAP  = 2.394632D+00
C         -1S-
      C( 1)= 8.93226887D-05
      C( 2)=-3.42831198D-04
      C( 3)= 6.00478236D-04
      C( 4)=-7.20350098D-04
      C( 5)= 2.54142307D-03
      C( 6)= 4.80099255D-02
      C( 7)= 2.97193922D-01
      C( 8)= 3.94011442D-01
      C( 9)= 2.41476753D-01
      C(10)= 1.06854677D-01
      C(11)= 4.02974468D-02
      C(12)= 1.44107767D-02
      C(13)= 4.89453604D-03
      C(14)= 1.73415135D-03
      C(15)= 5.27476375D-04
      C(16)= 2.37928933D-04
      C(17)= 2.64854313D-05
      C(18)= 4.87751540D-05
C         -2S-
      C(19)=-1.96291731D-03
      C(20)= 9.07101465D-03
      C(21)= 1.56908586D-02
      C(22)= 4.32135503D-01
      C(23)= 6.08635070D-01
      C(24)= 9.35786270D-02
      C(25)=-1.91200714D-01
      C(26)=-1.65923366D-01
      C(27)=-8.04313271D-02
      C(28)=-3.23502848D-02
      C(29)=-1.16441050D-02
      C(30)=-4.10292663D-03
      C(31)=-1.37672969D-03
      C(32)=-4.88493710D-04
      C(33)=-1.47543587D-04
      C(34)=-6.68903520D-05
      C(35)=-7.37307712D-06
      C(36)=-1.36838651D-05
C         -2P-
      C(37)= 2.37317088D-04
      C(38)= 1.18479677D-03
      C(39)= 2.73326171D-03
      C(40)= 8.49019117D-02
      C(41)= 3.62235730D-01
      C(42)= 4.09698282D-01
      C(43)= 2.30455294D-01
      C(44)= 8.60242084D-02
      C(45)= 2.40356277D-02
      C(46)= 6.33594327D-03
      C(47)= 1.20355199D-03
      C(48)= 4.71710486D-04
C         -3S-
      C(49)= 3.46144280D-01
      C(50)= 6.61587046D-01
      C(51)= 2.55031296D-01
      C(52)=-3.73292448D-01
      C(53)=-2.81638356D-01
      C(54)=-3.22437539D-02
      C(55)= 6.26378695D-02
      C(56)= 5.24322976D-02
      C(57)= 2.46764271D-02
      C(58)= 9.88411610D-03
      C(59)= 3.52370888D-03
      C(60)= 1.24688463D-03
      C(61)= 4.15245708D-04
      C(62)= 1.48422544D-04
      C(63)= 4.44202464D-05
      C(64)= 2.03204362D-05
      C(65)= 2.20542921D-06
      C(66)= 4.14535367D-06
C         -3P-
      C(67)= 1.79139451D-01
      C(68)= 4.28354175D-01
      C(69)= 4.36849606D-01
      C(70)= 1.30257752D-01
      C(71)=-1.14113035D-01
      C(72)=-1.22279182D-01
      C(73)=-6.62306589D-02
      C(74)=-2.38722719D-02
      C(75)=-6.57930778D-03
      C(76)=-1.70870194D-03
      C(77)=-3.27026033D-04
      C(78)=-1.26648537D-04
      GO TO 5000
C
C           ARGON
C
  800 CONTINUE
      ALPHAS = 8.019604D-02
      BETAS  = 2.371994D+00
      ALPHAP = 4.900201D-02
      BETAP  = 2.385538D+00
C         -1S-
      C( 1)= 6.06523885D-05
      C( 2)=-2.50139707D-04
      C( 3)= 3.79185807D-04
      C( 4)=-3.49832854D-04
      C( 5)= 2.08013465D-03
      C( 6)= 5.58597537D-02
      C( 7)= 3.08205505D-01
      C( 8)= 3.89079373D-01
      C( 9)= 2.33740937D-01
      C(10)= 1.02880581D-01
      C(11)= 3.87726414D-02
      C(12)= 1.39131813D-02
      C(13)= 4.73747219D-03
      C(14)= 1.68888380D-03
      C(15)= 5.13523027D-04
      C(16)= 2.34651228D-04
      C(17)= 2.52766050D-05
      C(18)= 4.85819917D-05
C         -2S-
      C(19)=-1.94405164D-03
      C(20)= 9.53251124D-03
      C(21)= 1.90272358D-02
      C(22)= 4.47853140D-01
      C(23)= 6.05126445D-01
      C(24)= 7.34170797D-02
      C(25)=-1.97108985D-01
      C(26)=-1.63883614D-01
      C(27)=-7.85146459D-02
      C(28)=-3.14863463D-02
      C(29)=-1.13420517D-02
      C(30)=-4.01048244D-03
      C(31)=-1.34989796D-03
      C(32)=-4.81735354D-04
      C(33)=-1.45533498D-04
      C(34)=-6.68006200D-05
      C(35)=-7.13094458D-06
      C(36)=-1.38040027D-05
C         -2P-
      C(37)= 1.69894084D-05
      C(38)= 2.96888975D-04
      C(39)= 2.22282275D-03
      C(40)= 8.15532666D-02
      C(41)= 3.62114755D-01
      C(42)= 4.12311299D-01
      C(43)= 2.29608117D-01
      C(44)= 8.52513537D-02
      C(45)= 2.37628281D-02
      C(46)= 6.28530868D-03
      C(47)= 1.19162796D-03
      C(48)= 4.74347500D-04
C         -3S-
      C(49)= 3.42380324D-01
      C(50)= 6.62585556D-01
      C(51)= 2.66174332D-01
      C(52)=-3.89957738D-01
      C(53)=-2.90865200D-01
      C(54)=-2.55936098D-02
      C(55)= 6.66942608D-02
      C(56)= 5.35244553D-02
      C(57)= 2.48744976D-02
      C(58)= 9.94243988D-03
      C(59)= 3.54386609D-03
      C(60)= 1.25969194D-03
      C(61)= 4.20313454D-04
      C(62)= 1.51300424D-04
      C(63)= 4.52196035D-05
      C(64)= 2.09762598D-05
      C(65)= 2.19812876D-06
      C(66)= 4.32048110D-06
C         -3P-
      C(67)= 1.67741090D-01
      C(68)= 4.14816220D-01
      C(69)= 4.47256567D-01
      C(70)= 1.49264619D-01
      C(71)=-1.16419957D-01
      C(72)=-1.28758074D-01
      C(73)=-6.85695287D-02
      C(74)=-2.45924159D-02
      C(75)=-6.74053329D-03
      C(76)=-1.76114829D-03
      C(77)=-3.34850226D-04
      C(78)=-1.32269620D-04
C
C        GENERATE THE EXPONENTS
C
 5000 CONTINUE
      E(1) = ALPHAS*BETAS
      DO I=2,NS
         E(I) = E(I-1) * BETAS
      ENDDO
      E(NS+1) = ALPHAP*BETAP
      DO I=2,NP
         E(I+NS) = E(I+NS-1) * BETAP
      ENDDO
C
C        IN CASE WE ARE FLOATING THE OUTER GAUSSIANS
C
      IF(NUCZ.LE.12) THEN
         C(64) = 1.0D+00
         RETURN
      ELSE
         C(79) = 1.0D+00
         C(80) = 1.0D+00
      END IF
      RETURN
      END
C
C*MODULE VVOS    *DECK BNDANX
      SUBROUTINE BNDANX
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXAO=8192, MXATM=2000, MXRT=100)
      DOUBLE PRECISION LOCMO
      DIMENSION FANT(8),LFANT(8),GANT(27),LGANT(8),LOCMO(3),
     *          NRNFG(10),NPFLG(10)
      LOGICAL GOPARR,DSKWRK,MASWRK,OK,ABELSYM,ABEL,ABELPT,
     *        ANALYS,FDIRCT,QCORR,DIRSCF,FDIFF,DIRTRF
C
      COMMON /DETPAR/ ICLBBR,ANALYS
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,S,SZ,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORB,
     *                NADET,NBDET,K,KST,IROOT,IPURES,MAXW1,NITER,
     *                MAXP,NCI,IGPDET,KSTSYM,NFTGCI
      COMMON /FCCWFN/ NSPACE,MSTA(51),MNUM(51),MINI(51),MAXI(51),
     *                IAMI(51),IAMA(51),IBMI(51),IBMA(51),IDIM(51),
     *                LBST(51),NREF0,FDIRCT,QCORR,C0SQ
      COMMON /FMCOM / XX(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NEX,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /PRPOPT/ ILOCAL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMMOL/ GROUP,COMPLEX,IGROUP,NAXIS,ILABMO,ABELSYM
      COMMON /VVOPAR/ IVVOS,NVVOS,BNDDEN
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C 
C      ----- INITIALIZE $BNDANA INPUT GROUP ----- 
C 
      PARAMETER (NNAM=1)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA BNDGRP/8HBNDANA  /
      DATA   QNAM/8HDENS    /
      DATA  KQNAM/5/
C
      DATA SCF,CISD,CISDT,CISDTQ,ORMAS
     *   /8HSCF     ,8HCISD    ,8HCISDT   ,8HCISDTQ  ,8HORMAS   /
      DATA LOCMO/8HBOYS    ,8HRUEDNBRG,8HPOP     /
      DATA CHECK/8HCHECK   /
C
      DATA FANT/8HC1      ,8HCI      ,8HCS      ,8HC2      ,
     *          8HD2      ,8HC2V     ,8HC2H     ,8HD2H     /
      DATA LFANT/1,1,1,1,2,2,2,3/
      DATA GANT/8HA       ,8HAG      ,8HAU      ,8HA'      ,
     *          8HA"      ,8HA       ,8HB       ,8HA       ,
     *          8HB1      ,8HB2      ,8HB3      ,8HA1      ,
     *          8HA2      ,8HB1      ,8HB2      ,8HAG      ,
     *          8HBG      ,8HBU      ,8HAU      ,8HAG      ,
     *          8HB1G     ,8HB2G     ,8HB3G     ,8HAU      ,
     *          8HB1U     ,8HB2U     ,8HB3U     /
      DATA LGANT/0,1,3,5,7,11,15,19/
C
      DATA ONE,TWO/1.0D+00,2.0D+00/
      DATA NRNFG /1,1,0,1,1,0,0,0,0,0/
      DATA NPFLG /0,0,0,0,0,0,0,0,0,0/
C
      DMATRX=CISD
      JRET = 0
      CALL NAMEIO(IR,JRET,BNDGRP,NNAM,QNAM,KQNAM,
     *            DMATRX,  0,0,0,
     *   0,0,0,0,0,    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *   0,0,0,0,0,    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *   0,0,0,0,0,    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0)
      IF(JRET.EQ.2) THEN
         IF(MASWRK) WRITE(IW,9000)
         CALL ABRT
      END IF
C
      OK = .FALSE.
      IF(DMATRX.EQ.SCF)    OK=.TRUE.
      IF(DMATRX.EQ.CISD)   OK=.TRUE.
      IF(DMATRX.EQ.CISDT)  OK=.TRUE.
      IF(DMATRX.EQ.CISDTQ) OK=.TRUE.
      IF(DMATRX.EQ.ORMAS)  THEN
         IF(CITYP.EQ.ORMAS) THEN
            OK=.TRUE.
         ELSE
            IF(MASWRK) WRITE(IW,9010)
         END IF
      END IF
C
      IF(OK) THEN
         BNDDEN = DMATRX
      ELSE
         IF(MASWRK) WRITE(IW,9020) DMATRX,SCF,CISD,CISDT,CISDTQ,ORMAS
         CALL ABRT
      END IF
C
C        --- PREPARE A CLOSED SHELL WAVEFUNCTION, AND ITS VVOS ---
C        IN CASE CITYP=ORMAS IS USED, THE NEXT CALL ALSO GENERATES
C        THE DENSITY MATRIX TO BE USED.
C
      CALL ENERGX
C
      MCORE = NUMCOR()
      IF(MASWRK) WRITE(IW,9030) LOCMO(ILOCAL),SCFTYP,
     *                          NA-MCORE,NVVOS,DMATRX
C
C       --- GET A DENSITY FOR THESE SEMICANONICAL VALENCE ORBITALS ---
C
      IF(DMATRX.EQ.SCF) THEN
         IF(MASWRK) WRITE(IW,9040)
      END IF
C
      IF(DMATRX.EQ.ORMAS) THEN
         IF(MASWRK) WRITE(IW,9050)
      END IF
C
C           USUAL OPTION IS A HARDWIRED GROUND STATE CI-SD WITH ORMAS
C           IT IS TRIVIAL TO EXTEND THIS TO HIGHER EXCITATIONS.
C
      IF(DMATRX.EQ.CISD    .OR.  
     *   DMATRX.EQ.CISDT   .OR.
     *   DMATRX.EQ.CISDTQ) THEN
C
         IF(DMATRX.EQ.CISD)   IEXCIT=2
         IF(DMATRX.EQ.CISDT)  IEXCIT=3
         IF(DMATRX.EQ.CISDTQ) IEXCIT=4
C
         SZ     = (MUL-1)/TWO
         IPURES = 1
C
         NCOR   = MCORE
         NCORSV = MCORE
         NACT   = NA - MCORE + NVVOS
         NORB   = MCORE + NACT
C
         NHIGH = INT(SZ+SZ+0.0001D+00)
         NELS  = NA + NB - 2*MCORE
         NBDET = (NELS-NHIGH)/2
         NADET = NBDET+NHIGH
C
         KSTSYM=1
         GRPDET=FANT(1)
         IF(IGROUP.EQ.1)                GRPDET = FANT(1)
         IF(IGROUP.EQ.3)                GRPDET = FANT(2)      
         IF(IGROUP.EQ.2)                GRPDET = FANT(3)
         IF(IGROUP.EQ.4.AND.NAXIS.EQ.2) GRPDET = FANT(4)
         IF(IGROUP.EQ.8.AND.NAXIS.EQ.2) GRPDET = FANT(5)      
         IF(IGROUP.EQ.7.AND.NAXIS.EQ.2) GRPDET = FANT(6)      
         IF(IGROUP.EQ.6.AND.NAXIS.EQ.2) GRPDET = FANT(7)      
         IF(IGROUP.EQ.9.AND.NAXIS.EQ.2) GRPDET = FANT(8)
         DO I=1,8
            IF (GRPDET.EQ.FANT(I)) THEN
               IGPDET=LFANT(I)
               STSYM = GANT(LGANT(I)+KSTSYM)
            ENDIF
         END DO
C
         IROOT=1
         CALL VCLR(WSTATE,1,MXRT)
         WSTATE(1) = ONE
         CALL VICLR(NFLGDM,1,MXRT)
         NFLGDM(1) = 1
         CALL VICLR(IWTS,1,MXRT)
         IWTS(1)=1
C
         MAXW1  = 300    ! NHGSS
         CRIT   = 1D-06  ! CVGTOL
         NITER  = 100    ! ITERMX
         K      = 1      ! NSTATE
         KST    = 1      ! NSTGSS
         MAXP   = 10     ! MXXPAN
         PRTTOL = 0.05D+00
C
         ICLBBR=0
         ANALYS=.FALSE.
C
         IF(MASWRK) WRITE(IW,9060) IEXCIT,NCORSV,NACT,NA-NCORSV,
     *                             NVVOS,NELS
C
C           DEFINE TWO ORMAS SPACES: NA (OF INFOA) INCLUDES CORE E-
C
         NSPACE  = 2
         MSTA(1) = NCORSV+1
         MSTA(2) = NA+1
         MSTA(3) = NA+NVVOS+1
         MNUM(1) = NA-NCORSV
         MNUM(2) = NVVOS
         MINI(1) = NELS-IEXCIT
         MINI(2) = 0
         MAXI(1) = NELS
         MAXI(2) = IEXCIT
         FDIRCT  = .TRUE.
         QCORR   = .FALSE.
         CALL FCCHECK(IW,.FALSE.,
     *                NSPACE,MNUM,MINI,MAXI,IAMI,IAMA,IBMI,IBMA,
     *                NADET,NBDET)
C
                    DIRTRF=DIRSCF
         IF(GOPARR) DIRTRF=.TRUE.
C
         ABEL = ABELPT()
         IF(.NOT.ABEL  .AND.  .NOT.DIRTRF) THEN
            IF(MASWRK) WRITE(IW,9070)
            CALL SYMOFF
            CALL JANDK
            CALL SYMON
         END IF
C
C           PREPARE ORBITAL SYMMETRY LABELS
C
         L0 = NQMT
         L1 = NUM
         L2 = (L1*L1+L1)/2
         L3 =  L1*L1
C
         CALL VALFM(LOADFM)
         LMOLAB = LOADFM + 1
         LMOIRP = LMOLAB + L1
         LVEC   = LMOIRP + L1
         LS     = LVEC   + L3
         LQ     = LS     + L2
         LWRK   = LQ     + L3
         LMODEG = LWRK   + L1
         LAST   = LMODEG + L1
         NEEDL  = LAST - LOADFM - 1
         CALL GETFM(NEEDL)
C
         IF(GRPDET.EQ.FANT(1)  .OR.  EXETYP.EQ.CHECK) THEN
            CALL C1DET(XX(LMOIRP),XX(LMOLAB),L0)
         ELSE
            CALL DAREAD(IDAF,IODA,XX(LVEC),L3,15,0)
            CALL DAREAD(IDAF,IODA,XX(LS),L2,12,0)
            CALL DAREAD(IDAF,IODA,XX(LQ),L3,45,0)
            CALL TRFSYM(XX(LMOLAB),XX(LMOIRP),XX(LMODEG),XX(LQ),XX(LS),
     *                  XX(LVEC),XX(LWRK),IA,L0,L1,L0,L1)
         END IF
         CALL GAJASW(XX(LMOIRP),NUM,GRPDET)
         CALL DAWRIT(IDAF,IODA,XX(LMOIRP),L1,262,1)
         CALL RETFM(NEEDL)
C
C           CARRY OUT THE CI-SD/SDT/SDTQ COMPUTATION WITH ORMAS.
C
         CALL ORDET(NRNFG,NPFLG)
      END IF
C
C        --- OBTAIN THE LOCALIZATION TRANSFORMATION ONTO ATOMS ---
C
      CALL LMOX
C
C        --- ROTATE DENSITY MATRIX TO THE LOCALIZED ORBITAL BASIS ---
C        -M1- INCLUDES OCCUPIED AND VIRTUAL VALENCE, BUT NOT CORE
C
      MDOC = NA - MCORE
      M1 = MDOC + NVVOS
      M2 = (M1*M1+M1)/2
      M3 = M1*M1
C
      CALL VALFM(LOADFM)
      LDCMO = LOADFM + 1
      LDLMO = LDCMO  + M2
      LTRAN = LDLMO  + M2
      LWRK  = LTRAN  + M3
      LAST  = LWRK   + M1
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
      CALL BNDANADEN(DMATRX,XX(LDCMO),XX(LDLMO),XX(LTRAN),XX(LWRK),
     *               M1,M2,MDOC)
      CALL RETFM(NEED)
C
C        --- ORIENT THE QUASIATOMIC ORBITALS AND PRINT ANALYSIS ---
C
      CALL DIRLMO
C
      RETURN
C
 9000 FORMAT(1X,'ERROR READING $BNDANA GROUP, TRY AGAIN')
 9010 FORMAT(1X,'*** ERROR ***'/
     *       1X,'THE USE OF DENS=ORMAS REQUIRES CITYP=ORMAS, AS WELL'/
     *       1X,'AS $CIDET AND $ORMAS INPUT GROUPS TO DEFINE THE CI.')
 9020 FORMAT(1X,'INVALID CHOICE FOR DENSITY MATRIX=',A8/
     *       1X,'PLEASE CHOOSE ONLY FROM: ',A8,4(',',A8),'.')
 9030 FORMAT(/5X,24(1H=)/5X,'BONDING ANALYSIS DETAILS'/5X,24(1H=)/
     *     1X,'THE ANALYSIS WILL BE BASED ON ',A8,
     *        ' LOCALIZATION OF ',A8,' ORBITALS,'/
     *     1X,'WHICH COMPRISE',I5,' OCCUPIED VALENCE AND',I5,
     *        ' VIRTUAL VALENCE ORBITALS,'/
     *     1X,'USING A DENSITY MATRIX OBTAINED AT THE ',A6,' LEVEL.')
 9040 FORMAT(/1X,'USING THE SIMPLISTIC SCF LEVEL DENSITY IS NOT AS',
     *           ' REVEALING AS A CI DENSITY.')
 9050 FORMAT(/1X,'THE CUSTOMIZED CI COMPUTATION IS ALREADY FINISHED,',
     *           ' AFTER THE SCF CONVERGED.')
 9060 FORMAT(/1X,'PERFORMING THE CI CALCULATION AT EXCITATION',
     *           ' LEVEL',I2,', USING ORMAS,'/
     *        1X,'WITH',I5,' CHEMICAL CORE AND',I5,' VALENCE ORBITALS',
     *           ' (OCC=',I5,', VIRT=',I5,'),'/
     *        1X,'CORRELATING ALL',I5,' VALENCE ELECTRONS.')
 9070 FORMAT(/1X,'REGENERATING AO INTEGRALS WITHOUT EXPLOITING POINT',
     *           ' GROUP SYMMETRY'/
     *        1X,'FOR USE BY THE NON-ABELIAN GROUP INTEGRAL',
     *           ' TRANSFORMATION')
      END
C
C*MODULE VVOS    *DECK BNDANADEN
      SUBROUTINE BNDANADEN(DMATRX,DCMO,DLMO,TRAN,WRK,M1,M2,MDOC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION DCMO(M2),DLMO(M2),TRAN(M1,M1),WRK(M1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA TWO/2.0D+00/
      DATA SCF/8HSCF     /
      DATA CHECK/8HCHECK   /
C
C         PREPARE THE USE OF JOE'S ORIENTATION PROGRAM BY GETTING
C         THE CORRECT DENSITY MATRIX TO THE CORRECT DISK RECORD.
C
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C         READ THE TRANSFORMATION TO LOCALIZED ORBITALS
C
      CALL DAREAD(IDAF,IODA,TRAN,M1*M1,73,0)
C
C         READ (OR GENERATE) THE DENSITY MATRIX
C
      IF(DMATRX.EQ.SCF) THEN
         CALL VCLR(DCMO,1,M2)
         DO I=1,MDOC
            II = (I*I+I)/2
            DCMO(II) = TWO
         ENDDO
      ELSE
         CALL DAREAD(IDAF,IODA,DCMO,M2,320,0)
      END IF
C
C          ROTATE THE DENSITY TO THE LOCALIZED ORBITAL BASIS
C
      CALL TFTRI(DLMO,DCMO,TRAN,WRK,M1,M1,M1)
C
C          WRITE OUT THE LOCALIZED DENSITY TO DAF RECORD 285
C
  800 CONTINUE
      IF(EXETYP.EQ.CHECK) CALL VCLR(DLMO,1,M2)
      CALL DAWRIT(IDAF,IODA,DLMO,M2,285,0)
      RETURN
      END
