C 08 AUG 07 - JHS - ADD SUBROUTINE NUCIDEN TO GET INITIAL DENSITY FOR PROPERTY RUNS
C 25 OCT 06 - SPW - NEW MODULE FOR CALCULATION OF NEO-CI DENSITY MATRICES
C
C*MODULE NEODEN  *DECK NEODEN_ELNU 
      SUBROUTINE NEODEN_ELNU(IECOUP,INCOUP,CICOEF,DMXNEO,
     *                DM1,DM2,DM1_NU,DM2_NU,INDEX,INDEXE,INDEXN,
     *                IA1,IB1,IA2,IB2,NALP,NBLP,NACT,NA,NB,
     *                IA1_NU,IB1_NU,IA2_NU,IB2_NU,NALP_NU,NBLP_NU,
     *                NACT_NU,NA_NU,NB_NU,
     *                M1,M2,M1_NU,M2_NU,IA,MX2,NCI,ISTATE,
     *                DOTWON)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DOTWON
C
      DIMENSION IECOUP(*),INCOUP(*),CICOEF(*),DMXNEO(*)
      DIMENSION DM1(*),DM2(*),DM1_NU(*),DM2_NU(*) 
      DIMENSION INDEX(M2,M2_NU),INDEXE(M2,M2),INDEXN(M2_NU,M2_NU),IA(*)
      DIMENSION IA1(NA),IB1(NA),IA2(NA),IB2(NA)
      DIMENSION IA1_NU(NA_NU),IB1_NU(NA_NU),IA2_NU(NA_NU),IB2_NU(NA_NU)
C
C     --- SUBROUTINE TO RETURN THE NEO-CI MIXED EL-NU DENSITY MATRIX
C         AS WELL AS THE 1-PARTICLE AND 2-PARTICLE DENSITIES ---
C
C                   WRITTEN BY SIMON P. WEBB          
C
C     ON ENTRY: CICOEF  - NUCLEAR-ELECTRONIC CI COEFICIENTS
C               IECOUP  - 1-E COUPLING COEFFICENTS
C               INCOUP  - 1-N COUPLING COEFFICENTS
C               INDEX   - ELECTRONIC-NUCLEAR INTEGRAL INDEXING
C               INDEXE  - ELECTRONIC INTEGRAL INDEXING
C               INDEXN  - NUCLEAR INTEGRAL INDEXING
C               ISTATE  - CI ROOT # 
C
C     ON EXIT:  DMXNEO - NUCLEAR-ELECTRONIC CI DENSITY MATRIX  
C               DM1    - ELECTRONIC  1-PARTICLE CI DENSITY MATRIX  
C               DM2    - ELECTRONIC  2-PARTICLE CI DENSITY MATRIX  
C               DM1_NU - NUCLEAR  1-PARTICLE CI DENSITY MATRIX  
C               DM2_NU - NUCLEAR  2-PARTICLE CI DENSITY MATRIX  
C
C
C     --- SET UP INDEXING ARRAY (IJ=ACTIVE NUC; KL=ACTIVE ELECS) ---
C
      DO KL=1,M2
         DO  IJ=1,M2_NU
            INDEX(KL,IJ) = (IJ-1)*M2 + KL
         END DO
      END DO
C
      DO IJ=1,M2
         DO  KL=1,IJ
            INDEXE(IJ,KL) = IJ*(IJ-1)/2 + KL
            INDEXE(KL,IJ) = INDEXE(IJ,KL)
         END DO
      END DO
C
      DO IJ=1,M2_NU
         DO  KL=1,IJ
            INDEXN(IJ,KL) = IJ*(IJ-1)/2 + KL
            INDEXN(KL,IJ) = INDEXN(IJ,KL)
         END DO
      END DO
C
      DO II=1,MX2
         IA(II) = (II*II+II)/2
      END DO
C
      IHELEM = 0
C
      DO 20 I=1,NA_NU
         IA1_NU(I) = I
   20 CONTINUE
C
C     --- OUTER LOOP OVER ALPHA NUCLEAR DETERMINANTS ---
C
      IMIX1 = 1
      I1_NU = 0
      DO 30 IIA=1,NALP_NU
C
         DO 40 I=1,NB_NU
            IB1_NU(I) = I
   40    CONTINUE
C
C        --- OUTER LOOP OVER BETA NUCLEAR DETERMINANTS ---
C
         DO 50 IIB=1,NBLP_NU
            I1_NU = I1_NU + 1
C
            I1 = 0
            DO 60 I=1,NA
               IA1(I) = I
   60       CONTINUE
C
C           --- OUTER LOOP OVER ALPHA ELECTRONIC DETERMINANTS ---
C
            DO 70 KKA=1,NALP
C
               DO 80 I=1,NB
                  IB1(I) = I
   80          CONTINUE
C
C              --- OUTER LOOP OVER BETA ELECTRONIC DETERMINANTS ---
C
               DO 100 KKB=1,NBLP
                  I1 = I1 + 1
C
                  I2_NU = 0
                  DO 110 I=1,NA_NU
                     IA2_NU(I) = I
  110             CONTINUE
C
C                 --- INNER LOOP OVER ALPHA NUCLEAR DETERMINANTS ---
C
                  IMIX2 = 1
                  DO 120 JJA=1,NALP_NU
C
                     DO 130 I=1,NB_NU
                        IB2_NU(I) = I
  130                CONTINUE
C  
C                    --- LOOP OVER BETA NUCLEAR DETERMINANTS ---
C
                     DO 140 JJB=1,NBLP_NU
                        I2_NU = I2_NU + 1
C
                        I2 = 0
                        DO 150 I=1,NA
                           IA2(I) = I
  150                   CONTINUE
C  
C                       --- INNER LOOP OVER ALPHA ELECTRONIC DETERMINANTS ---
C
                        DO 160 LLA=1,NALP
C
                           DO 170 I=1,NB
                              IB2(I) = I
  170                      CONTINUE
C
C                          -- INNER LOOP OVER BETA ELECTRONIC DETERMINANTS --
C
                           DO 180 LLB=1,NBLP
                              IF(IMIX2.GT.IMIX1) GOTO 350
                              I2 = I2 + 1
C
C                             --- GET NUC-ELEC CROSS TERM DENSITY ELEM ---
C
                              IHELEM = IHELEM + 1
                              CALL NUELDN(DMXNEO,CICOEF,INDEX,
     *                             NA,NB,IA1,IB1,I1,I2,
     *                             NA_NU,NB_NU,IA1_NU,IB1_NU,
     *                             I1_NU,I2_NU,
     *                             M2,M2_NU,IA,IECOUP,INCOUP,NCI,
     *                             IMIX1,IMIX2,ISTATE)
C
C                             --- GET ELEC CI DENSITY MATRIX ELEMENTS ---
C
                              IF(I1_NU.EQ.I2_NU) THEN
                                 CALL DEN12_NEO(DM1,DM2,CICOEF,M1,
     *                                   NA,NB,IA1,IB1,
     *                                   IA2,IB2,0,I1,I2,INDEXE,
     *                                   IMIX1,IMIX2,
     *                                   NCI,ISTATE)
                              END IF
C
C                             --- GET NUC CI DENSITY MATRIX ELEMENTS ---
C
                              IF(I1.EQ.I2 .AND. DOTWON) THEN
                                 CALL DEN12_NEO(DM1_NU,DM2_NU,CICOEF,
     *                                M1_NU,NA_NU,NB_NU,IA1_NU,IB1_NU,
     *                                IA2_NU,IB2_NU,0,I1_NU,I2_NU,
     *                                INDEXN,
     *                                IMIX1,IMIX2,NCI,ISTATE)
                              ELSE IF(I1.EQ.I2) THEN
                                 CALL DEN1_NEO(DM1_NU,CICOEF,
     *                                M1_NU,NA_NU,NB_NU,IA1_NU,IB1_NU,
     *                                IA2_NU,IB2_NU,0,I1_NU,I2_NU,
     *                                INDEXN,
     *                                IMIX1,IMIX2,NCI,ISTATE)
                              END IF
C
C                             --- LOOP BACK FOR NEXT ELEMENTS --- 
C
                              CALL ADVANC(IB2,NB,NACT)
                              IMIX2 = IMIX2 + 1
  180                      CONTINUE
C
                           CALL ADVANC(IA2,NA,NACT)
C
  160                   CONTINUE
C
                        CALL ADVANC(IB2_NU,NB_NU,NACT_NU)
C
  140                CONTINUE
C
                     CALL ADVANC(IA2_NU,NA_NU,NACT_NU)
C
  120             CONTINUE
C
  350             CONTINUE
                  CALL ADVANC(IB1,NB,NACT)
                  IMIX1 = IMIX1 + 1
C
  100          CONTINUE
C
               CALL ADVANC(IA1,NA,NACT)
C
   70       CONTINUE
C
            CALL ADVANC(IB1_NU,NB_NU,NACT_NU)
C
   50    CONTINUE
         CALL ADVANC(IA1_NU,NA_NU,NACT_NU)
C
   30 CONTINUE
C
      RETURN
      END
C
C*MODULE NEODEN  *DECK NUELDN
      SUBROUTINE NUELDN(DMXNEO,CICOEF,INDEX,
     *                  NA,NB,IA1,IB1,I1,I2,
     *                  NA_NU,NB_NU,IA1_NU,IB1_NU,
     *                  I1_NU,I2_NU,M2,M2_NU,IDXT,IECOUP,
     *                  INCOUP,NCI,IMIX1,IMIX2,NRT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DMXNEO(*),CICOEF(NCI,*),INDEX(M2,M2_NU),IDXT(*)
      DIMENSION IA1(NA),IB1(NB)
C     DIMENSION IA2(NA),IB2(NB)
      DIMENSION IA1_NU(NA_NU),IB1_NU(NB_NU)
C     DIMENSION IA2_NU(NA_NU),IB2_NU(NB_NU)
      DIMENSION IECOUP(*),INCOUP(*)
C
C     RETURNS THE MATRIX ELEMENT < IK | H(MIX) |J L > WHERE K, L ARE
C     DETERMINANTS.  ALPHA AND BETA OCCUPATIONS ARE STORED IN
C     IACON1, IBCON1 FOR K AND IACON2, IBCON2 FOR L.
C
      IF(I1.GT.I2) THEN
         IE  = I1*(I1-1)/2 + I2
      ELSE
         IE  = I2*(I2-1)/2 + I1
      END IF 
      IF(I1_NU.GT.I2_NU) THEN
         IN = I1_NU*(I1_NU-1)/2 + I2_NU
      ELSE
         IN = I2_NU*(I2_NU-1)/2 + I1_NU
      END IF 
C
C     --- IF ELEC DETS THE SAME AND NUC DETS THE SAME ---
C
      IF (I1.EQ.I2 .AND. I1_NU.EQ.I2_NU) THEN
C
C        --- ALPHA ELECS - ALPHA NUC ---
         DO 10 I=1,NA
            IA = IDXT(IA1(I))
            DO 20 J=1,NA_NU
              JA = IDXT(IA1_NU(J))
              IX = INDEX(IA,JA)
              DMXNEO(IX)=DMXNEO(IX)+CICOEF(IMIX1,NRT)*CICOEF(IMIX2,NRT)
   20       CONTINUE
   10    CONTINUE
C
C        --- BETA ELECS - ALPHA NUC ---
         DO 30 I=1,NB
            IB = IDXT(IB1(I))
            DO 40 J=1,NA_NU
              JA = IDXT(IA1_NU(J))
              IX = INDEX(IB,JA)
              DMXNEO(IX)=DMXNEO(IX)+CICOEF(IMIX1,NRT)*CICOEF(IMIX2,NRT) 
   40       CONTINUE
   30    CONTINUE
C
C        --- ALPHA ELECS - BETA NUC ---
         DO 50 I=1,NA
            IA = IDXT(IA1(I))
            DO 60 J=1,NB_NU
              JB = IDXT(IB1_NU(J))
              IX = INDEX(IA,JB)
              DMXNEO(IX)=DMXNEO(IX)+CICOEF(IMIX1,NRT)*CICOEF(IMIX2,NRT) 
   60       CONTINUE
   50    CONTINUE
C
C        --- BETA ELECS - BETA NUC ---
         DO 70 I=1,NB
            IB = IDXT(IB1(I))
            DO 80 J=1,NB_NU
              JB = IDXT(IB1_NU(J))
              IX = INDEX(IB,JB)
              DMXNEO(IX)=DMXNEO(IX)+CICOEF(IMIX1,NRT)*CICOEF(IMIX2,NRT) 
   80       CONTINUE
   70    CONTINUE
C
C
C     --- IF ELECTRONIC DETERMINANTS ARE SAME ---
C
      ELSE IF (I1.EQ.I2) THEN
C
         NUIDX = INCOUP(IN)
         IF(NUIDX.EQ.0) GOTO 120
         JA = ABS(NUIDX)
         SIGN = DBLE(NUIDX/JA)
C
C        --- ALPHA ELECS ---
         DO 100 I=1,NA
            IA = IDXT(IA1(I))
            IX = INDEX(IA,JA)
            DMXNEO(IX) = DMXNEO(IX)
     *                 + SIGN*CICOEF(IMIX1,NRT)*CICOEF(IMIX2,NRT)
  100    CONTINUE
C
C        --- BETA ELECS ---
         DO 110 I=1,NB
            IB = IDXT(IB1(I))
            IX = INDEX(IB,JA)
            DMXNEO(IX) = DMXNEO(IX) 
     *                 + SIGN*CICOEF(IMIX1,NRT)*CICOEF(IMIX2,NRT)
  110    CONTINUE
C
  120    CONTINUE
C
C     --- IF NUCLEAR DETERMINANTS ARE SAME ---
C
      ELSE IF (I1_NU.EQ.I2_NU) THEN
C
         IDXEL = IECOUP(IE)
         IF(IDXEL.EQ.0) GOTO 160
         IA = ABS(IDXEL)
         SIGN = DBLE(IDXEL/IA)
C
C        --- ALPHA NUC ---
         DO 140 J=1,NA_NU
            JA = IDXT(IA1_NU(J))
            IX = INDEX(IA,JA)
            DMXNEO(IX) = DMXNEO(IX)
     *                 + SIGN*CICOEF(IMIX1,NRT)*CICOEF(IMIX2,NRT)
  140    CONTINUE
C
C        ---  BETA NUC ---
         DO 150 J=1,NB_NU
            JB = IDXT(IB1_NU(J))
            IX = INDEX(IA,JB)
            DMXNEO(IX) = DMXNEO(IX)
     *                 + SIGN*CICOEF(IMIX1,NRT)*CICOEF(IMIX2,NRT)
  150    CONTINUE
C
  160    CONTINUE
C
C     --- BOTH ELEC AND NUCLEAR DETERMINANTS DIFFERENT ---
C
      ELSE 
C
         NUIDX = INCOUP(IN)
         IF(NUIDX.EQ.0) GOTO 200
         JA = ABS(NUIDX)
         SIGN1 = DBLE(NUIDX/JA)
C
         IDXEL = IECOUP(IE)
         IF(IDXEL.EQ.0) GOTO 200
         IA = ABS(IDXEL)
         SIGN2 = DBLE(IDXEL/IA)
C
         IX = INDEX(IA,JA)
         DMXNEO(IX) = DMXNEO(IX)
     *              + SIGN1*SIGN2*CICOEF(IMIX1,NRT)*CICOEF(IMIX2,NRT)
C
  200    CONTINUE
C
      ENDIF
C
      RETURN
      END
C
C*MODULE NEODEN  *DECK DEN12_NEO
      SUBROUTINE DEN12_NEO(DM1,DM2,CICOEF,NORB,NALPHA,NBETA,IA1,IB1,
     *                     IA2,IB2,NCOR,IJ,JI,INDEX,
     *                     IMIX1,IMIX2,NCI,IIST)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION IA1(NALPHA),IB1(NBETA)
      DIMENSION IA2(NALPHA),IB2(NBETA)
      DIMENSION DM1(*),DM2(*),CICOEF(NCI,*)
      DIMENSION INDEX((NORB*(NORB+1))/2,*)
      INTEGER DIFF1(2),DIFF2(2),IPOS1(2),IPOS2(2)
C
C     ALPHA AND BETA OCCUPATIONS ARE STORED IN
C     IA1, IB1 FOR K AND IA2, IB2 FOR L.
C
      DIFF1(1) = 0
      DIFF1(2) = 0
      DIFF2(1) = 0
      DIFF2(2) = 0
C
C     --- IF DETERMINANTS ARE THE SAME ---
C
      IF (IJ.EQ.JI) THEN
C
         DO 200 I=1,NALPHA
            IA = IA1(I)
            I1 = INDEX(IA,IA)
            DM1(I1) = DM1(I1) + CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
C
            DO 197 J=I+1,NALPHA
               IIA1 = IA1(J)
               I2 = INDEX(IIA1,IIA1)
               IT = INDEX(I1,I2)
               DM2(IT) = DM2(IT) + CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
               IC = INDEX(IA,IIA1)
               IT = INDEX(IC,IC)
               DM2(IT) = DM2(IT) - CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
  197       CONTINUE
C
            DO 198 J=1,NBETA
               IIB1 = IB1(J)
               I2 = INDEX(IIB1,IIB1)
               IT = INDEX(I1,I2)
               DM2(IT) = DM2(IT) + CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
  198       CONTINUE
  200    CONTINUE
C
         DO 210 I=1,NBETA
            IB = IB1(I)
            I1 = INDEX(IB,IB)
            DM1(I1) = DM1(I1) + CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
            DO 204 J=I+1,NBETA
               IIB1 = IB1(J)
               I2 = INDEX(IIB1,IIB1)
               IT = INDEX(I1,I2)
               DM2(IT) = DM2(IT) + CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST) 
               IC = INDEX(IB,IIB1)
               IT = INDEX(IC,IC)
               DM2(IT) = DM2(IT) - CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
  204       CONTINUE
  210    CONTINUE
         RETURN
C
      END IF
C
      IDEA=0
      IDEB=0
C
C     --- DIFFERENT ORBITALS IN FIRST DETERMINANT ---
C
      DO 20 I=NCOR+1,NALPHA
         DO 15 J=NCOR+1,NALPHA
            IF (IA1(I).EQ.IA2(J)) GOTO 20
   15    CONTINUE
         IDEA = IDEA + 1
         IF (IDEA.GT.2) RETURN
         DIFF1(IDEA) = IA1(I)
         IPOS1(IDEA) = I
   20 CONTINUE
C
      DO 30 I=NCOR+1,NBETA
         DO 25 J=NCOR+1,NBETA
            IF (IB1(I).EQ.IB2(J)) GOTO 30
   25    CONTINUE
         IDEB = IDEB + 1
         IF (IDEA+IDEB.GT.2) RETURN
         DIFF1(IDEA+IDEB) = IB1(I)
         IPOS1(IDEA+IDEB) = I
   30 CONTINUE
C
C     --- FIND THE DIFFERENT ORBITALS IN SECOND DETERMINANAT ---
C
      IST = NCOR+1
      DO 63 II=1,IDEA
            DO 50 I=IST,NALPHA
               DO 45 J=NCOR+1,NALPHA
                  IF (IA2(I).EQ.IA1(J)) GOTO 50
   45          CONTINUE
               GOTO 60
   50       CONTINUE
C
   60       DIFF2(II) = IA2(I)
            IPOS2(II) = I
            IST = I+1
   63 CONTINUE
C
      IST = NCOR+1
      DO 163 II=1,IDEB
            DO 150 I=IST,NBETA
               DO 145 J=NCOR+1,NBETA
                  IF (IB2(I).EQ.IB1(J)) GOTO 150
  145          CONTINUE
               GOTO 160
  150       CONTINUE
  160       DIFF2(II+IDEA) = IB2(I)
            IPOS2(II+IDEA) = I
            IST = I+1
  163 CONTINUE
C
C     --- IF DETERMINANTS DIFFER BY ONE ORBITAL ---
C
      IF (IDEA+IDEB.EQ.1) THEN
C
C        --- ONE PARTICLE DENSITY CONTRIBUTION ---
C
         IND1 = INDEX(DIFF1(1),DIFF2(1))
         IPERM = (-1)**(IPOS1(1)-IPOS2(1))
         DM1(IND1) = DM1(IND1) +
     *               IPERM*CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
C
C        --- TWO PARTICLE DENSITY CONTRIBUTION ---
C
C        --- IF DIFFERENT ORBITALS ARE ALPHA SPIN ORBS ---
C
         IF (IDEA.EQ.1) THEN
            DO 673 K=1,NALPHA
               NK = IA1(K)
               IF (NK.EQ.DIFF1(1)) GOTO 673
               IND2 = INDEX(NK,NK)
               INDX = INDEX(IND1,IND2)
               DM2(INDX) = DM2(INDX) +
     *                     CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)*IPERM
               I1 = INDEX(DIFF1(1),NK)
               I2 = INDEX(DIFF2(1),NK)
               INX = INDEX(I1,I2)
               DM2(INX) = DM2(INX) -
     *                    IPERM*CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
  673       CONTINUE
C
            DO 678 K=1,NBETA
               NK = IB1(K)
               IND2 = INDEX(NK,NK)
               INDX = INDEX(IND1,IND2)
               DM2(INDX) = DM2(INDX) +
     *                     IPERM*CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
  678       CONTINUE
C
C        --- DIFFERENT ORBITALS ARE BETA SPIN ORBS ---
C
         ELSE
C
            DO 732 K=1,NALPHA
               NK = IA1(K)
               IND2 = INDEX(NK,NK)
               INDX = INDEX(IND1,IND2)
               DM2(INDX) = DM2(INDX) +
     *                     IPERM*CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
  732       CONTINUE
C
            DO 752 K=1,NBETA
               NK = IB1(K)
               IF (NK.EQ.DIFF1(1)) GOTO 752
               IND2 = INDEX(NK,NK)
               INDX = INDEX(IND1,IND2)
               DM2(INDX) = DM2(INDX) +
     *                     IPERM*CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
               I1 = INDEX(DIFF1(1),NK)
               I2 = INDEX(DIFF2(1),NK)
               INX = INDEX(I1,I2)
               DM2(INX) = DM2(INX) -
     *                    IPERM*CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
  752       CONTINUE
C
         END IF
         RETURN
C
      ELSE
C
C        TWO ORBITALS ARE DIFFERENT
C        CONTRIBUTION ONLY TO 2-PARTICLE DENSITY MATRIX.
C        DIFFERING ORBITALS IN DIFF1(1),DIFF1(2) FOR CON1 AND
C        DIFF2(1),DIFF2(2) FOR CON2.  POSITION STORED IN
C        IPOS1(1),IPOS1(2) AND IPOS2(1),IPOS2(2).
C
         IPERM = (-1)**(IPOS1(1)-IPOS2(1)+IPOS1(2)-IPOS2(2))
         I11 = DIFF1(1)
         I12 = DIFF2(1)
         I21 = DIFF1(2)
         I22 = DIFF2(2)
         CALL FLSHBF(6)
         I1 = INDEX(I11,I12)
         I2 = INDEX(I21,I22)
         INX = INDEX(I1,I2)
         DM2(INX) = DM2(INX) +
     *              IPERM*CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
C
C        --- IF ALL DIFFERING ORBITALS ARE OF SAME SPIN THEN ---
C            HAVE EXTRA MATRIX ELEMENTS.
C
         IF (IDEA.EQ.2.OR.IDEB.EQ.2) THEN
            I1 = INDEX(I11,I22)
            I2 = INDEX(I12,I21)
            INX = INDEX(I1,I2)
            DM2(INX) = DM2(INX) -
     *                 IPERM*CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
         ENDIF
C
      ENDIF
C
      RETURN
      END
C
C*MODULE NEODEN  *DECK NEONO
      SUBROUTINE NEONO(DMO,VAO,VNO,OCCNO,IWRK,WRK,SCR,
     *                 NCORE,NACT,M1,M2,L1,COROCC,IW)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DMO(M2),VAO(L1,L1),VNO(L1,L1),OCCNO(L1),
     *          IWRK(M1),WRK(M1,8),SCR(M1)
C
      PARAMETER (ZERO=0.0D+00)
C
C     --- GENERATE NATURAL ORBITALS ---
C
      NOCC = NCORE + NACT
      NVIRT = L1-NOCC
C
C     --- DIAGONALIZE MO DENSITY MATRIX ---
C
      CALL GLDIAG(L1,M1,M1,DMO,WRK,SCR,VNO(1,NCORE+1),IERR,IWRK)
      IF(IERR.NE.0) THEN
         WRITE(IW,*) 'DIAG FAILED IN DETNO'
         CALL ABRT
      END IF
C
C     --- REORDER FROM BIGGEST TO LOWEST EIGENVALUES ---
C
      NDUM = M1/2
      DO 140 IDUM = 1,NDUM
         NFIRST = IDUM
         NLAST = M1-IDUM+1
         DUM = SCR(NLAST)
         SCR(NLAST ) = SCR(NFIRST)
         SCR(NFIRST) = DUM
         CALL DSWAP(M1,VNO(1,NCORE+NFIRST),1,VNO(1,NCORE+NLAST),1)
  140 CONTINUE
C
C     --- TRANSFORM ACTIVE NO'S TO ATOMIC BASIS ---
C
      CALL TFSQB(VNO(1,NCORE+1),VAO(1,NCORE+1),OCCNO,NACT,L1,L1)
C
C     --- COPY ANY CORE AND VIRTUAL ORBITALS, ADJUST PHASE OF NOS ---
C
      IF(NCORE.GT.0) CALL DCOPY(L1*NCORE,VAO,1,VNO,1)
      IF(NVIRT.GT.0) CALL DCOPY(L1*NVIRT,VAO(1,NOCC+1),1,
     *                                   VNO(1,NOCC+1),1)
      CALL STFASE(VNO,L1,L1,L1)
C
C     --- SET UP OCCUPATION NUMBER ARRAY ---
C
      DO 210 I = 1,NCORE
         OCCNO(I) = COROCC
  210 CONTINUE
      DO 220 I = 1,NACT
         OCCNO(I+NCORE) = SCR(I)
  220 CONTINUE
      DO 230 I=1,NVIRT
         OCCNO(I+NOCC) = ZERO
  230 CONTINUE
C
      RETURN
C
      END
C
C*MODULE NEODEN  *DECK NEODM1
      SUBROUTINE NEODM1(CICOEF,
     *                DM1,DM1_NU,INDEXE,INDEXN,
     *                IA1,IB1,IA2,IB2,NALP,NBLP,NACT,NA,NB,
     *                IA1_NU,IB1_NU,IA2_NU,IB2_NU,NALP_NU,NBLP_NU,
     *                NACT_NU,NA_NU,NB_NU,
     *                M1,M2,M1_NU,M2_NU,IA,MX2,NCI,NROOT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION CICOEF(*)
      DIMENSION DM1(*),DM1_NU(*)
      DIMENSION INDEXE(M2,M2),INDEXN(M2_NU,M2_NU),IA(*)
      DIMENSION IA1(NA),IB1(NA),IA2(NA),IB2(NA)
      DIMENSION IA1_NU(NA_NU),IB1_NU(NA_NU),IA2_NU(NA_NU),IB2_NU(NA_NU)
C
C     --- SUBROUTINE TO RETURN THE NEO-CI ELECTRONIC AND NUCLEAR 
C                                ONE-PARTICLE DENSITY MATRICES ONLY --- 
C
C                    WRITTEN BY SIMON P. WEBB
C
C     ON ENTRY: CICOEF  - NUCLEAR-ELECTRONIC CI COEFICIENTS
C               IECOUP  - 1-E COUPLING COEFFICENTS
C               INCOUP  - 1-N COUPLING COEFFICENTS
C               INDEXE  - ELECTRONIC INTEGRAL INDEXING
C               INDEXN  - NUCLEAR INTEGRAL INDEXING
C               NROOT   - CI ROOT # 
C
C     ON EXIT:  DM1    - ELECTRONIC  1-PARTICLE CI DENSITY MATRIX  
C               DM1_NU - NUCLEAR  1-PARTICLE CI DENSITY MATRIX  
C
C
C     --- SET UP INDEXING ARRAY (IJ=ACTIVE NUC; KL=ACTIVE ELECS) ---
C
      DO IJ=1,M2
         DO  KL=1,IJ
            INDEXE(IJ,KL) = IJ*(IJ-1)/2 + KL
            INDEXE(KL,IJ) = INDEXE(IJ,KL)
         END DO
      END DO
C
      DO IJ=1,M2_NU
         DO  KL=1,IJ
            INDEXN(IJ,KL) = IJ*(IJ-1)/2 + KL
            INDEXN(KL,IJ) = INDEXN(IJ,KL)
         END DO
      END DO
C
      DO II=1,MX2
         IA(II) = (II*II+II)/2
      END DO
C
      DO 20 I=1,NA_NU
         IA1_NU(I) = I
   20 CONTINUE
C
C     --- OUTER LOOP OVER ALPHA NUCLEAR DETERMINANTS ---
C
      IMIX1 = 1
      I1_NU = 0
      DO 30 IIA=1,NALP_NU
C
         DO 40 I=1,NB_NU
            IB1_NU(I) = I
   40    CONTINUE
C
C        --- OUTER LOOP OVER BETA NUCLEAR DETERMINANTS ---
C
         DO 50 IIB=1,NBLP_NU
            I1_NU = I1_NU + 1
C
            I1 = 0
            DO 60 I=1,NA
               IA1(I) = I
   60       CONTINUE
C
C           --- OUTER LOOP OVER ALPHA ELECTRONIC DETERMINANTS ---
C
            DO 70 KKA=1,NALP
C
               DO 80 I=1,NB
                  IB1(I) = I
   80          CONTINUE
C
C              --- OUTER LOOP OVER BETA ELECTRONIC DETERMINANTS ---
C
               DO 100 KKB=1,NBLP
                  I1 = I1 + 1
C
                  I2_NU = 0
                  DO 110 I=1,NA_NU
                     IA2_NU(I) = I
  110             CONTINUE
C
C                 --- INNER LOOP OVER ALPHA NUCLEAR DETERMINANTS ---
C
                  IMIX2 = 1
                  DO 120 JJA=1,NALP_NU
C
                     DO 130 I=1,NB_NU
                        IB2_NU(I) = I
  130                CONTINUE
C  
C                    --- LOOP OVER BETA NUCLEAR DETERMINANTS ---
C
                     DO 140 JJB=1,NBLP_NU
                        I2_NU = I2_NU + 1
C
                        I2 = 0
                        DO 150 I=1,NA
                           IA2(I) = I
  150                   CONTINUE
C  
C                       --- INNER LOOP OVER ALPHA ELECTRONIC DETERMINANTS ---
C
                        DO 160 LLA=1,NALP
C
                           DO 170 I=1,NB
                              IB2(I) = I
  170                      CONTINUE
C
C                          -- INNER LOOP OVER BETA ELECTRONIC DETERMINANTS --
C
                           DO 180 LLB=1,NBLP
                              IF(IMIX2.GT.IMIX1) GOTO 350
                              I2 = I2 + 1
C
C                             --- GET ELEC CI DENSITY MATRIX ELEMENTS ---
C
                              IF(I1_NU.EQ.I2_NU) THEN
                                 CALL DEN1_NEO(DM1,CICOEF,M1,
     *                                   NA,NB,IA1,IB1,
     *                                   IA2,IB2,0,I1,I2,INDEXE,
     *                                   IMIX1,IMIX2,
     *                                   NCI,NROOT)
                              END IF
C
C                             --- GET NUC CI DENSITY MATRIX ELEMENTS ---
C
                              IF(I1.EQ.I2) THEN
                                 CALL DEN1_NEO(DM1_NU,CICOEF,
     *                                M1_NU,NA_NU,NB_NU,IA1_NU,IB1_NU,
     *                                IA2_NU,IB2_NU,0,I1_NU,I2_NU,
     *                                INDEXN,
     *                                IMIX1,IMIX2,NCI,NROOT)
                              END IF
C
C                             --- LOOP BACK FOR NEXT ELEMENTS --- 
C
                              CALL ADVANC(IB2,NB,NACT)
                              IMIX2 = IMIX2 + 1
  180                      CONTINUE
C
                           CALL ADVANC(IA2,NA,NACT)
C
  160                   CONTINUE
C
                        CALL ADVANC(IB2_NU,NB_NU,NACT_NU)
C
  140                CONTINUE
C
                     CALL ADVANC(IA2_NU,NA_NU,NACT_NU)
C
  120             CONTINUE
C
  350             CONTINUE
                  CALL ADVANC(IB1,NB,NACT)
                  IMIX1 = IMIX1 + 1
C
  100          CONTINUE
C
               CALL ADVANC(IA1,NA,NACT)
C
   70       CONTINUE
C
            CALL ADVANC(IB1_NU,NB_NU,NACT_NU)
C
   50    CONTINUE
         CALL ADVANC(IA1_NU,NA_NU,NACT_NU)
C
   30 CONTINUE
C
      RETURN
      END
C
C*MODULE NEODEN  *DECK DEN1_NEO
      SUBROUTINE DEN1_NEO(DM1,CICOEF,NORB,NALPHA,NBETA,IA1,IB1,
     *                     IA2,IB2,NCOR,IJ,JI,INDEX,
     *                     IMIX1,IMIX2,NCI,IIST)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION IA1(NALPHA),IB1(NBETA)
      DIMENSION IA2(NALPHA),IB2(NBETA)
      DIMENSION DM1(*),CICOEF(NCI,*)
      DIMENSION INDEX((NORB*(NORB+1))/2,*)
      INTEGER DIFF1(2),DIFF2(2),IPOS1(2),IPOS2(2)
C
C     ALPHA AND BETA OCCUPATIONS ARE STORED IN
C     IA1, IB1 FOR K AND IA2, IB2 FOR L.
C
      DIFF1(1) = 0
      DIFF1(2) = 0
      DIFF2(1) = 0
      DIFF2(2) = 0
C
C     --- IF DETERMINANTS ARE THE SAME ---
C
      IF (IJ.EQ.JI) THEN
C
         DO 200 I=1,NALPHA
            IA = IA1(I)
            I1 = INDEX(IA,IA)
            DM1(I1) = DM1(I1) + CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
  200    CONTINUE
C
         DO 210 I=1,NBETA
            IB = IB1(I)
            I1 = INDEX(IB,IB)
            DM1(I1) = DM1(I1) + CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
  210    CONTINUE
         RETURN
C
      END IF
C
      IDEA=0
      IDEB=0
C
C     --- DIFFERENT ORBITALS IN FIRST DETERMINANT ---
C
      DO 20 I=NCOR+1,NALPHA
         DO 15 J=NCOR+1,NALPHA
            IF (IA1(I).EQ.IA2(J)) GOTO 20
   15    CONTINUE
         IDEA = IDEA + 1
         IF (IDEA.GT.2) RETURN
         DIFF1(IDEA) = IA1(I)
         IPOS1(IDEA) = I
   20 CONTINUE
C
      DO 30 I=NCOR+1,NBETA
         DO 25 J=NCOR+1,NBETA
            IF (IB1(I).EQ.IB2(J)) GOTO 30
   25    CONTINUE
         IDEB = IDEB + 1
         IF (IDEA+IDEB.GT.2) RETURN
         DIFF1(IDEA+IDEB) = IB1(I)
         IPOS1(IDEA+IDEB) = I
   30 CONTINUE
C
C     --- FIND THE DIFFERENT ORBITALS IN SECOND DETERMINANAT ---
C
      IST = NCOR+1
      DO 63 II=1,IDEA
            DO 50 I=IST,NALPHA
               DO 45 J=NCOR+1,NALPHA
                  IF (IA2(I).EQ.IA1(J)) GOTO 50
   45          CONTINUE
               GOTO 60
   50       CONTINUE
C
   60       DIFF2(II) = IA2(I)
            IPOS2(II) = I
            IST = I+1
   63 CONTINUE
C
      IST = NCOR+1
      DO 163 II=1,IDEB
            DO 150 I=IST,NBETA
               DO 145 J=NCOR+1,NBETA
                  IF (IB2(I).EQ.IB1(J)) GOTO 150
  145          CONTINUE
               GOTO 160
  150       CONTINUE
  160       DIFF2(II+IDEA) = IB2(I)
            IPOS2(II+IDEA) = I
            IST = I+1
  163 CONTINUE
C
C     --- IF DETERMINANTS DIFFER BY ONE ORBITAL ---
C
      IF (IDEA+IDEB.EQ.1) THEN
C
C        --- ONE PARTICLE DENSITY CONTRIBUTION ---
C
         IND1 = INDEX(DIFF1(1),DIFF2(1))
         IPERM = (-1)**(IPOS1(1)-IPOS2(1))
         DM1(IND1) = DM1(IND1) +
     *               IPERM*CICOEF(IMIX1,IIST)*CICOEF(IMIX2,IIST)
C
      ELSE
C
C     --- IF DETERMINANTS DIFFER BY MORE NO CONTRIBUTION --- 
C
         RETURN
C
      ENDIF
C
      RETURN
      END
C
C*MODULE NEODEN  *DECK NEODM2
      SUBROUTINE NEODM2(IECOUP,INCOUP,CICOEF,DMXNEO,
     *                DM1,DM2,DM1_NU,DM2_NU,INDEX,INDEXE,INDEXN,
     *                IA1,IB1,IA2,IB2,NALP,NBLP,NACT,NA,NB,
     *                IA1_NU,IB1_NU,IA2_NU,IB2_NU,NALP_NU,NBLP_NU,
     *                NACT_NU,NA_NU,NB_NU,
     *                M1,M2,M4,M1_NU,M2_NU,M4_NU,IA,MX2,MIXDET,
     *                NSTATE,WSTATE,SADMX,SADM1,SADM2,
     *                SADM1_NU,SADM2_NU,DOTWON)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DOTWON
C
      DIMENSION IECOUP(*),INCOUP(*),CICOEF(*),DMXNEO(*)
      DIMENSION DM1(*),DM2(*),DM1_NU(*),DM2_NU(*) 
      DIMENSION INDEX(M2,M2_NU),INDEXE(M2,M2),INDEXN(M2_NU,M2_NU),IA(*)
      DIMENSION IA1(NA),IB1(NA),IA2(NA),IB2(NA)
      DIMENSION IA1_NU(NA_NU),IB1_NU(NA_NU),IA2_NU(NA_NU),IB2_NU(NA_NU)
      DIMENSION WSTATE(NSTATE)
      DIMENSION SADMX(*),SADM1(*),SADM2(*),SADM1_NU(*),SADM2_NU(*) 
C
C     --- DRIVER FOR THE NEO-CI MIXED EL-NU DENSITY MATRIX
C         AS WELL AS THE 1-PARTICLE AND 2-PARTICLE DENSITIES
C         AND STATE-AVERAGING ---
C
      CALL VCLR(SADMX,1,M2*M2_NU)
      CALL VCLR(SADM1,1,M2)
      CALL VCLR(SADM2,1,M4)
      CALL VCLR(SADM1_NU,1,M2_NU)
      CALL VCLR(SADM2_NU,1,M4_NU)
C
      DO 100 ISTATE = 1,NSTATE
         WEI = WSTATE(ISTATE)
         IF(WEI.LT.1.0D-12) GO TO 100
         CALL VCLR(DMXNEO,1,M2*M2_NU)
         CALL VCLR(DM1,1,M2)
         CALL VCLR(DM2,1,M4)
         CALL VCLR(DM1_NU,1,M2_NU)
         CALL VCLR(DM2_NU,1,M4_NU)
C
         ICIRT = ISTATE
C
         CALL NEODEN_ELNU(IECOUP,INCOUP,CICOEF,DMXNEO,
     *        DM1,DM2,DM1_NU,DM2_NU,INDEX,INDEXE,INDEXN,
     *        IA1,IB1,IA2,IB2,NALP,NBLP,NACT,NA,NB,
     *        IA1_NU,IB1_NU,IA2_NU,IB2_NU,NALP_NU,NBLP_NU,
     *        NACT_NU,NA_NU,NB_NU,
     *        M1,M2,M1_NU,M2_NU,IA,MX2,MIXDET,
     *        ICIRT,DOTWON)
C
C        --- STATE-AVERAGED ELECTRONIC DENSITY MATRICES ---
C
         DO I=1,M2
            SADM1(I) = SADM1(I) +  WEI*DM1(I)
         END DO
C
         DO I=1,M4
            SADM2(I) = SADM2(I) +  WEI*DM2(I)
         END DO
C
C        --- STATE-AVERAGED NUCLEAR DENSITY MATRICES ---
C
         DO I=1,M2_NU
            SADM1_NU(I) = SADM1_NU(I) +  WEI*DM1_NU(I)
         END DO
C
         DO I=1,M4_NU
            SADM2_NU(I) = SADM2_NU(I) +  WEI*DM2_NU(I)
         END DO
C
C        --- STATE-AVERAGED MIXED DENSITY MATRIX ---
C
         DO I=1,M2*M2_NU
            SADMX(I)  = SADMX(I)  +  WEI*DMXNEO(I)
         END DO
C
  100 CONTINUE
C
      RETURN
      END
C
C*MODULE NEODEN  *DECK WTDNEO
      SUBROUTINE WTDNEO(EXETYP,DM1,DM2,LBABEL,X,IX,NINTMX,LABSIZ,
     *                  M1,M2,M4,WRK,NOCC2,NCORE,CUTOFF,
     *                  NFTO,NRECO,NDM2O)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL IEQJ,KEQL,GOPARR,DSKWRK,MASWRK
C
      DIMENSION DM1(M2),DM2(M4),LBABEL(M1),
     *          X(NINTMX),IX(*),WRK(NOCC2)
      DIMENSION MULT8(8),LKUPSM(64)
C
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /IOFILE/ IR,IW,IP,IS,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (HALF=0.5D+00, TWO=2.0D+00)
C
      DATA MULT8/0,8,16,24,32,40,48,56/
      DATA LKUPSM/1,2,3,4,5,6,7,8,
     *            2,1,4,3,6,5,8,7,
     *            3,4,1,2,7,8,5,6,
     *            4,3,2,1,8,7,6,5,
     *            5,6,7,8,1,2,3,4,
     *            6,5,8,7,2,1,4,3,
     *            7,8,5,6,3,4,1,2,
     *            8,7,6,5,4,3,2,1/
      DATA CHECK/8HCHECK   /
C
C     ----- WRITE -DM1- AND -DM2- -----
C
      THRSH = 1.0D+01*CUTOFF
C
      DO 50 I=1,M1
         IF(LBABEL(I).EQ.0) THEN
            IBAD = I+NCORE
            IF(MASWRK) WRITE(IW,9000) IBAD
            CALL ABRT
         END IF
   50 CONTINUE
C
      CALL DSCAL(M4,HALF*HALF,DM2,1)
C
      SMALL = 1.0D-07
      IF(EXETYP.EQ.CHECK) SMALL = -HALF
      NERR = 0
      IJ = 0
      DO 110 I=1,M1
         ISYM = LBABEL(I)
         DO 100 J=1,I
            IJ = IJ+1
            IJMUL = MULT8(ISYM)+LBABEL(J)
            IJSYM = LKUPSM(IJMUL)
            IF(IJSYM.NE.1) THEN
               IF(ABS(DM1(IJ)).LT.THRSH) THEN
                  DM1(IJ) = 0.0D+00
               ELSE
                  IBAD = I+NCORE
                  JBAD = J+NCORE
                  IF(MASWRK) WRITE(IW,9010) IBAD,JBAD,DM1(IJ)
                  IF(MASWRK) WRITE(IW,9030)
                  CALL ABRT
               END IF
            END IF
  100    CONTINUE
         IF(DM1(IJ).LT.SMALL) THEN
            IF(MASWRK) WRITE(IW,9040) I,DM1(IJ)
            NERR=NERR+1
         END IF
  110 CONTINUE
C
      IF(NERR.GT.0) THEN
         IF(MASWRK) WRITE(IW,*) 'CHECK YOUR ACTIVE SPACE CAREFULLY.'
         IF(MASWRK) WRITE(IW,*) 'THE 1ST ORDER DENSITY MATRIX IS:'
         CALL PRTRI(DM1,M1)
         CALL ABRT
      END IF
C
C     ----- WRITE NONZERO ELEMENTS OF -DM2- ON FILE -NFTO- -----
C
      NRECO=0
      NDM2O=0
C
      N = 0
      NX = 1
      DO 280 I = 1,M1
         ISYM = LBABEL(I)
         DO 260 J = 1,I
            IEQJ = I.EQ.J
            IJMUL = MULT8(ISYM)+LBABEL(J)
            IJSYM = LKUPSM(IJMUL)
            DO 240 K = 1,I
               LMAX = K
               IF(K.EQ.I) LMAX = J
               IJKMUL = MULT8(IJSYM)+LBABEL(K)
               IJKSYM = LKUPSM(IJKMUL)
               DO 220 L = 1,LMAX
                  KEQL=K.EQ.L
                  N = N+1
                  VAL = DM2(N)
C
                  IF(IEQJ)                VAL=VAL+VAL
                  IF(KEQL)                VAL=VAL+VAL
                  IF(I.EQ.K .AND. J.EQ.L) VAL=VAL+VAL
                  IF(ABS(VAL).LT.CUTOFF) GO TO 220
C
C      ONLY TOTALLY SYMMETRIC DIRECT PRODUCT SHOULD BE NONZERO ELEMENTS
C
                  LSYM = LBABEL(L)
                  IF(LSYM.NE.IJKSYM) THEN
                     IF(ABS(VAL).LT.THRSH) GO TO 220
                     IBAD = I+NCORE
                     JBAD = J+NCORE
                     KBAD = K+NCORE
                     LBAD = L+NCORE
                     IF(MASWRK) WRITE(IW,9020) IBAD,JBAD,KBAD,LBAD,VAL
                     IF(MASWRK) WRITE(IW,9030)
                     CALL ABRT
                  END IF
C
              NPACK = NX
              IPACK = I
              JPACK = J
              KPACK = K
              LPACK = L
              IF (LABSIZ .EQ. 2) THEN
*I32            LABEL1 = ISHFT( IPACK, 16 ) + JPACK
*I32            LABEL2 = ISHFT( KPACK, 16 ) + LPACK
*I32            IX( 2*NPACK-1 ) = LABEL1
*I32            IX( 2*NPACK   ) = LABEL2
*I64            LABEL = ISHFT( IPACK, 48 ) + ISHFT( JPACK, 32 ) +
*I64 *                  ISHFT( KPACK, 16 ) + LPACK
*I64            IX(NPACK) = LABEL
              ELSE IF (LABSIZ .EQ. 1) THEN
*I32            LABEL = ISHFT( IPACK, 24 ) + ISHFT( JPACK, 16 ) +
*I32 *                  ISHFT( KPACK,  8 ) + LPACK
*I32            IX(NPACK) = LABEL
*I64            IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64              LABEL = ISHFT( IPACK, 24 ) + ISHFT( JPACK, 16 ) +
*I64 *                    ISHFT( KPACK,  8 ) + LPACK
*I64              IX( NPACK/2 ) = IX( NPACK/2 ) + LABEL
*I64            ELSE
*I64              LABEL = ISHFT( IPACK, 56 ) + ISHFT( JPACK, 48 ) +
*I64 *                    ISHFT( KPACK, 40 ) + ISHFT( LPACK, 32 )
*I64              IX( (NPACK/2)+1 ) = LABEL
*I64            END IF
              END IF
C
                  X(NX) = VAL
                  NX = NX+1
                  IF(NX.GT.NINTMX) THEN
                     CALL PWRIT(NFTO,X,IX,NINTMX,NINTMX)
                     NRECO=NRECO+1
                     NDM2O=NDM2O+NINTMX
                     NX = 1
                  END IF
  220          CONTINUE
  240       CONTINUE
  260    CONTINUE
  280 CONTINUE
C
C         OUTPUT THE FINAL, PARTIAL BUFFER OF DM2
C
      NX = -(NX-1)
      CALL PWRIT(NFTO,X,IX,NX,NINTMX)
      NDM2O=NDM2O-NX
      NRECO=NRECO+1
C
C     ----- OUTPUT FIRST ORDER DENSITY IN MO BASIS -----
C     GENERATE MO DENSITY OVER ALL ORBITALS, INCLUDING CORE
C
      CALL VCLR(WRK,1,NOCC2)
      II = 0
      DO 300 I=1,NCORE
         II = II+I
         WRK(II) = TWO
  300 CONTINUE
      DO 320 I=1,M1
         IV = IA(I)
         IC = IA(I+NCORE)
         DO 310 J=1,I
             IJV = IV + J
             IJC = IC + J + NCORE
             WRK(IJC) = DM1(IJV)
  310    CONTINUE
  320 CONTINUE
C
C     WRITE -DM1- WITHOUT CORE ORBITALS TO DAF RECORD 320
C     WRITE -DM1- WITH    CORE ORBITALS TO DAF RECORD 68
C     WRITE -DM1- WITH    CORE ORBITALS AFTER -DM2- ON FILE -NFTO-
C
      CALL DAWRIT(IDAF,IODA,DM1,   M2,320,0)
      CALL DAWRIT(IDAF,IODA,WRK,NOCC2, 68,0)
      CALL SQWRIT(NFTO,WRK,NOCC2)
      RETURN
C
 9000 FORMAT(1X,'UNABLE TO SIFT DENSITY MATRIX, ORBITAL',I5,
     *          ' HAS UNKNOWN SYMMETRY.')
 9010 FORMAT(/1X,'INACCURATE 1ST ORDER DENSITY MATRIX ELEMENT FOUND,'/
     *       1X,'GAMMA(',2I5,')=',E20.10,' FOUND, IT SHOULD BE ZERO'/
     *       1X,'BY SYMMETRY.')
 9020 FORMAT(/1X,'INACCURATE 2ND ORDER DENSITY MATRIX ELEMENT FOUND,'/
     *       1X,'GAMMA(',4I5,')=',E20.10,' FOUND, IT SHOULD BE ZERO'/
     *       1X,'BY SYMMETRY.')
 9030 FORMAT(/1X,'LOSS OF SYMMETRY IN THE DENSITY MATRIX MAY BE DUE TO'/
     *       1X,'    UNSYMMETRICAL ORBITALS: CHECK $VEC GROUP,'/
     *       1X,'         ADJUST $GUESS TOLZ=1.0D-5 TOLE=1.0D-04'/
     *       1X,'    OR, UNSYMMETRICAL CI ROOT: $DET CVGTOL=5.0D-06')
 9040 FORMAT(1X,'***** ERROR: ACTIVE ORBITAL',I3,
     *          ' HAS VERY SMALL OCCUPATION NUMBER=',1P,E13.6)
      END
C
C*MODULE NEO     *DECK NUCIDEN
      SUBROUTINE NUCIDEN
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      N1 = NUMNB
      N2 = (N1*N1+N1)/2
      N3 = N1*N1
C
C     --- GET NUCLEAR MO'S AND DENSITY ---
C
      CALL VALFM(LOADFM)
      LVA   = LOADFM + 1
      LAOC  = LVA    + N3
      LDA   = LAOC   + N1
      LAST  = LDA    + N2 
C
      NEED  = LAST   - LOADFM - 1
      CALL GETFM(NEED)
C
C
C     --- OCCUPATION NUMBERS ---
C
      DO 120 I = 1,N1
         X(I+LAOC-1) = ZERO
         IF(I.LE.NNA) X(I+LAOC-1) = ONE
  120 CONTINUE
C
C     --- GET NUCLEAR MO'S AND DENSITY ---
C
      CALL DAREAD(IDAF,IODA,X(LVA),N3,444,0)
      CALL DMTX(X(LDA),X(LVA),X(LAOC),NNA,N1,N1)
C
C     --- WRITE NUCLEAR DENSITY TO DAF ---
C
      CALL DAWRIT(IDAF,IODA,X(LDA),N2,445,0)
C
      CALL RETFM(NEED)
C
      RETURN
      END
